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DECISION BETWEEN TWO ALTERNATIVES--HOW MANY 
SXPERIMENTS?* 


P. M. Grunpy, D. HW. Rees M. J. R. Weary 


Rothamsted Experimental Station 


When in the course of scientific research a new technological process 
or agricultural treatment is suggested, the question at once arises 
whether the change-over from established practice is going to be 
economically justified. The answer to this question can be found by 
comparing the extra cost of the new process with the value of the extra 
output due to its adoption, and this latter quantity will have to be 
estimated by a set of experiments. Now experiments themselves cost 
money, and we are faced with a compromise; we can make certain of 
getting the right answer to our question at the expense of a very large 
experimental programme, or we can economise on our experiments at 
the risk of arriving at a wrong conclusion. It seems clear that there 
will be an optimum amount of experimentation which will in some 
sense give us the best return for our money, and this paper is concerned 
with the determination of such an optimum. 

The extra returns due to introducing the new process will depend 
on the seale on which it will be applied and on the increase in output 
per unit over the original process, and knowing these we could estimate 
the expected loss due to wrong decisions and balance it against the cost 
of experimentation. The increase in output per unit, however, is only 
determined from the results of the experiments themselves, and there 
are obvious advantages in some kind of sequential scheme in which the 
decision on whether or not to continue experimenting depends on the 
results so far obtained. With the present approach, the investigation 
of a fully sequential technique presents rather formidable difficulties, 
although some progress has been made in this direction; but a fairly 
detailed solution to the problem is possible under conditions analogous 
to double sampling as used in quality control. We suppose that a 
preliminary experiment is carried out, and that in the light of this we 
decide to adopt or reject the new process or else to postpone our decision 


*Read at the Third International Biometrie Congress, Bellagio, September, 1953. 


317 


pale 


318 BIOMETRICS, SEPTEMBER 1954 
until we have done one further set of experiments whose size depends 
on the results so far obtained. This is likely to be less efficient than a 
fully sequential scheme, but is of practical importance when a decision 
is needed quickly, so that there is only time to carry out two batches 
of experiments. A full mathematical treatment of the problem will be 
given in a forthcoming paper, the results of which are summarised here. 

We suppose, then, that the new process increases the output per 
unit by a quantity 7 and that the net gain due to introducing the new 
process* is given by 


k(n — 


where k’ is a known positive factor depending on the scale oi application 
and c is a known constant allowing for the difference in cost between 
the old and new processes and the capital cost of the change. The 
change-over will then be economically justifiable if (7 — c) is positive. 
To estimate 7 we carry out a set of experiments which provide a set of 
independent estimates y, normally distributed with standard deviation 
o. (We assume that o is known with sufficient accuracy from previous 
experience, and that it includes all sources of variation that affect the 
y’s.) Then, if we carry out n experiments, our decision as to the adoption 
or rejection of the new process will be based on the sign of (7 — c), 
where 7 is the mean of the n estimates y, and we can calculate the 
probability P of accepting the new process as a function of n. If the 
cost of each experiment is k, then the total cost of our experiments is 
kn and the expectation of gain (measured from the status quo) allowing 


for the uncertainty of our information is Py, so that the total of costs 
plus losses is 


R=kn — Pu 


which may be referred to as the total risk.The optimum value of n will 
be that which minimises R. However, we have still to face the difficulty 
that P, and hence R, is a function of the unknown 7. 

In the present approach, we carry out a preliminary experiment 
and obtain an estimate y,. P is taken to be the conditional probability 
of accepting the new process, given y, . If we now average R over the 
fiducial distribution of », which is normal with mean y, and standard 
deviation ¢, we arrive at an averaged risk R.** We choose the value of 


*If the gain is to continue over a period of years, a solution can be obtained by using an equivalent 
capital sum. 

**This approach seems intuitively reasonable, and results given below indicate that it provides a 
satisfactory solution to the present problem; but we do not wish to claim here any general optimal 
properties ‘for it. 
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n to minimise this quantity. On the present assumptions, this averaged 
risk takes the form 


R = kn — k'(y, — 


where 


#2) = [ dt. 


The possible modes of behaviour of F# as a function of n are illustrated 


diagramatically in Fig. 1. As n increases from 1, R may increase 


— > 


FIGURE 1. AVERAGED RISK AS A FUNCTION OF THE NUMBER OF EXPERIMENTS 


steadily or else may pass through a maximum and then a minimum 
value, thereafter increasing steadily. Furthermore, the minimum 
value may be greater than the value at n = 1. It is clear that only when 
the minimum value is less than the value at n = 1 is it worth while 
doing more experiments (the optimum number of experiments being 
given by the value of n at the minimum); otherwise, the decision to 
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adopt or reject the new process should be based on the result of the 
first experiment alone. 

The value of n which minimises R is actually given by the larger 
root of the equation 


This can be solved iteratively for given values of the constants; but it 
is much simpler to use the nomogram shown in Fig. 2. The left and 
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FIGURE 2. NOMOGRAM FOR DETERMINATION OF n, THE RECOMMENDED NUMBER 
OF EXPERIMENTS 


right hand scales are graduated in terms of 
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and a line joining the appropriate points on these scales intersects the 
curved scale at points corresponding to the two roots of the above 
equation. We have still to ensure that the minimum value of R is 
actually less than the value at n = 1. The condition for this to be true 
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ean be put in the simple form 


where 


and this has been incorporated in the nomogram in the form of the 
ungraduated curve AB. If the line joining the outer scales fails to 
intersect or to pass above this curve, the minimum R is greater than 
the value at n = 1 and no further experimentation is recommended. 
The recommended n can, alternatively, be obtained from Table 1. 


TABLE 1. RECOMMENDED VALUES OF n 


log Value of xr 
0 0.2 O 0.6 0.8 1.0 1.2 1.4 1.6 13s 20 23 24 

2.0 4.8 4.7 4.35 4.3 3.9 3.4 2.9 1 1 1 1 1 1 
2.2 5.9 5.8 5.7 5.3 4.9 4.3 3.7 3.0 1 1 1 1 1 
2.4 7.3 2. TF8 €7 6:3 5.5 4.8 4.0 1 1 1 1 1 
2.6 9.2 9.1 8.8 8.3 vat 7.0 6.1 5.2 4.1 1 1 1 1 
2.8 12 ll 11 10 9.7 8.8 7.7 6.6 D4 4.1 1 1 1 
3.0 14 14 14 13 12 11 9.7 8.4 7.0 5.35 4.0 1 1 
8.2 18 18 17 16 15 14 12 ll 8.9 7.2 5.5 1 1 
3.4 23 23 22 21 19 17 15 13 ll 9.3 7.3 5.3 | 
3.6 28 28 27 26 24 23 20 17 l4 12 9.5 7.2 1 


The values in this table were obtained from the nomogram, this me‘ hod 
being considered sufficiently accurate for the object in view. 
Having arrived at a rule for deciding on the amount of experimen- 
tation to be carried out, it remains to discuss the way in which the rule 
will operate in practice. It can be seen from Fig. 1 that, since the 
averaged risk always increases initially as n increases from 1 (except 
when 2 = 0 exactly), there is always a small amount of extra experi- 
mentation which is definitely uneconomic. This is very reasonable, 
since the only reason for doing more experiments is the possibility that 
they may contradict the initial result, and a fair amount of data will be 
necessary before we shall be prepared to go against our original findings. 
Fig. 3 shows the probability distribution of n, the optimum number of 
experiments, for \ = 10°°° and 6/¢ = 0, 1.5 and 3, where @ = 9 — c. 
Notice that whatever the result of the first experiment there is an upper 
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limit (depending on \) to the number of experiments that can be 
recommended. 


What the rule accomplishes can best be appreciated in terms of the 
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2 


FIGURE 3. FREQUENCY DISTRIBUTION OF n > 1 FOR 2X = 10%5 


expected cost plus losses under different circumstances. In Fig. 4, 
this total expected loss is plotted against @/o for \ equal to 107-° and 
10°. Also shown are the total losses for two “single sampling” schemes, 


where the amount of experimentation is decided once for all. The 
amounts considered are 
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where n,,. is the greatest number of experiments ever recommended 
under the double sampling rule. The second scheme gives n = 6.2 
ford = and n = 10.8 = 10°. The first scheme is uneconomical 
unless the effect of the new process is very small, when our potential 
gain to set against the cost of experimentation is small, or very large, 


Expected 
cost + 
loss 
ao} 
4 Azo? 
Expected 
cost + an 
loss 
Ral 
at 
aor 
set 
4.6 6/a 1 1 és. 
1 2 1 2 


FIGURE 4. EXPECTED (COSTS + LOSSES) FOR TWO SINGLE SAMPLING SCHEMES, 
n =1ANDn = in,,, AND FOR DOUBLE SAMPLING 


when a correct decision will nearly always be made; the second tends 
to be wasteful, especially when the effect of the new process turns out 
to be large. It is worth mentioning that no scheme on a singie or 
double sampling basis can give better results than the present rule for 
all values of 6/c. 

The original idea behind this research was given by Dr. F. Yates 
in an article in Nature (170, 138, 1952) where he gave a solution to the 


quantitative problem in which the optimum level of some treatm nt is 
to be determined. 
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THE ANALYSIS OF EXPERIMENTS CONTAINING 
DIFFERENT CROP ROTATIONS 


F. Yates 
Rothamsted Experimental Station. 


Summary 


The problems arising in the analysis of experiments containing 
different crop rotations are investigated. When the design of the 
experiment is such that each block coniains plots which sometimes 
carry a given crop but do not all carry the crop in the same set of years 
the year-block totals will not be orthogonal with the plot totals. In 
most such cases the fitting of constants must be resorted to in order 
to obtain separate estimates of plot error and plot X year error which 
are free of year X block interactions. The method is illustrated by 
application to a rice-pasture experiment containing rotations of different 
lengths and with different proportions of rice to pasture. 


Introduction 


The present paper deals only with problems arising in the analysis 
(not the design) of experiments containing different crop rotations, 
i.e. experiments of type (b) below, mainly in relation to a proposed 
experiment on alternative rice-pasture rotations in the United States. 
Some problems arising in the analysis of experiments comparing different 
treatments on the same rotation, i.e. experiments of type (a) below, 
have been discussed by Patterson (1953). A general discussion of the 
design of rotation experiments has been given by Yates (1949). An 
earlier discussion of the design and analysis of long-term experiments 
is provided by Cochran (1939), and some points arising in the analysis 


of experiments of type (b) are considered by Crowther and Cochran 
(1942). 


Terminology 


The design and analysis of experiments involving crop rotations is 
one that has received relatively little notice in the literature, and a 
brief note on the terminology of the subject may therefore be helpful 
to the understanding of this paper. 

Sequence. In this paper the term sequence is used to denote a 
sequence of crops, treatments or crop-treatment combinations which 
differs in any respect from other such sequences occurring in the ex- 
periment. 
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Cycle, period, phase. These terms are used in their customary sense. 
If we have a repetitive sequence of crops or treatments ¢, , ¢2 , ; 
Ca, Ci, Cz, ***, C, &@ Single repetition is termed a cycle. The number of 
years (or other time units) in the cycle is termed its period. Sequences 
starting at different points in the same cycle are said to be in different 
phase. There are thus n possible phases of a cycle of period n. The 
term phase is also used to denote the different components of the cycle. 
These alternative meanings are not likely to cause confusion and follow 
common usage (e.g. phases of an alternating electric current and phases 
of the moon). 

If two or more cycles of differing period are included in the experiment 
then the whole experiment will follow a cycle whose period is the lowest 
common multiple of the two or more periods. 


Rotation. A rotation is a definite cycle of crops grown in successive 
years on the same land. In agricultural practice the term is used 
somewhat loosely, and includes not only minor variations in cropping 
from cycle to cycle, e.g. the substitution of one cereal crop for another, 
but often also alterations in the length of the cycle. In experimental 
work such variations are naturally kept to a minimum. The separate 
crops of a rotation are in agricultural terminology called courses. If 
the same crop occurs twice in a rotation it will constitute two different 
courses. An n course rotation is therefore a cropping cycle of period n. 


Rotation experiments. There are two main classes of rotation ex- 
periments: 

(a) Experiments on the effects of treatments applied to a fixed 
rotation of crops. The experimental treatments may be repeated year 
after year, or may be varied in some manner which is rega: ed as 
appropriate to the questions at issue. A common device is t~ use a 
cycle of treatments. If the cycle of treatments on a given plot of a 
rotation experiment has the same period as the cropping cycle a given 
crop will have the same treatment each time it is grown. If the periods 
are different the treatment will vary in a cyclic manner, with the 
important consequence that any given combination of crop and treat- 
ment-phase occurs on different plots in successive crop cycles, thereby 


considerably increasing the accuracy. 


(b) Experiments comparing the effects of different rotations. Here 
the different crops themselves act as treatments, and plots between 
which comparisons have to be made will not always be carrying the 
same crop in the same year. This introduces the additional problem of 
design of ensuring that the plots to be compared simultaneously carry 
the same crop in a sufficient number of years for the necessary com- 
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parisons to be made, a problem which is particularly troublesome when 
the effect of rotations of different length is under investigation. 

Experiments of class (b) may, and often do, contain other experi- 
mental treatments in addition to the variations in cropping. 


Series. In rotation experiments of type (a) the area of land is usually 
divided into separate parts for the different crop phases. Except when 
the same crop occurs more than once in the cycle each part will in any 
one year carry a different crop. These separate parts are termed series. 

Blocks. The whole experiment, or if in series the separate series, 
will normally be divided into blocks, as in ordinary one-year agricultural 
experiments. Each block may contain a complete replicate of all the 
sequences in the experiment or the series, or the device of confounding 
may be used to reduce block size, in which case each block will only 
contain part of a replicate. 


Preliminary years. One, two or more years at the beginning of a 
cyclical experiment usually have to be excluded from the main analysis 
because treatments which would have been applied had the experiment 
been started earlier will not in fact have been applied. These years are 
known as the preliminary years. Their number will depend on the 
nature of the treatments and other factors. 


Phase differences. In rotation experiments many different types of 
contrast may be of interest. One type, which has no analogy in one- 
year experiments, is the contrast between different phases of the same 
cycle. Such contrasts may be termed phase differences. In a crop 
rotation containing the same crop twice in the rotation, for instance, 
there will be a phase difference for this crop. This phase difference can 
only be estimated with precision if the phases carrying the same crop 
occur in the same blocks. Thus with the crop rotation PQPR estima- 
tion of the phase difference requires arrangement in two series instead 
of four. With the rotation PPQR a single series must be used if estimates 
are to be available each year. 


The problem 


At the Summer Statistics Conference organized by the North 
Carolina Institute of Statistics in 1952 and held at Blue Ridge, North 
Carolina, an experiment to compare various rotations of rice and grass 
pasture was discussed. Comparisons were required between the 
rotations: 


A. 1 year rice, 2 years grass 
B. 2 years rice, 2 years grass 
C. 1 year rice, 3 years grass 
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The design proposed was to compare all phases of these three rotations 
in three randomised blocks of 11 plots each. A complete cycle of this 
experiment requires 12 years, after 3 preliminary years. One replicate 
of such a cycle (excluding the preliminary years) is shown in Table 1. 


TABLE 1. ONE REPLICATE OF THE RICE PASTURE ROTATION EXPERIMENT 


A letter denotes the plot is carrying rice. 


— 


Total 


1 a, bf Ci ¥; 
3 a3 bg obs C3 Y3 
4 a, bf by Y, 
5 as bs bf C5 Ys 
6 as bg be C6 Ye 
7 az bs b; C7 Y; 
8 ag 4 bs Cs 
9 ag by bg C9 Ys 
10 aio bfo bio Cio Yio 
ll an bf, bi Cu Yu 
12 ai2 bf. bie Ci2 Yi2 


Total P, P; P, P, Ps Py Pr Pu R 


The yields of rice in the three rotations can be assessed by comparing 
the means of the a, 3(b + b’) and c yields, and the difference between 
the first and second year’s rice in the B rotation (which we may ierm 
the crop phase difference) can be assessed by comparing the means of 
the b and b’ yields. All these comparisons will be free of year differences 
since each year is equally represented in each mean. The first group of 
comparisons will involve plot differences, whereas the phase comparisons 
will only involve plot and year interaction. Similar comparisons can 
be made over a shorter period than 12 years, though some of the 
symmetry is then lost. If only six years results are available, for 
example, each of the A plots is represented twice, but of the B plots 
plot 4 is represented four times, plots 5 and 7 three times, and plot 6 
twice, and similarly with the C plots. 

The problem is to estimate the errors of these and other comparisons 
that may require to be made. The ordinary subdivision of the analysis 
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of variance into a part derived from the totals of plots over all years, 
and another part (further sub-divided if necessary) derived from the 
plots X years interaction, breaks down, since the different plot totals 
involve year differences. 

The full analysis of an experiment of this type provides an interesting 
example of the partition of degrees of freedom in material which 
possesses a certain degree of balance and orderliness but is not fully 
orthogonal. In order to elucidate the various points at issue we will 
first consider the analysis of two somewhat simpler types of experiment. 


A rotation with all crops different, treatments tied to crops 


As a simple example we may take a three course rotation experiment 
to compare two treatment cycles (or fixed treatments), with three 
replicates of each of the three phases. Denote the treatment cycles by 
Aand B. It is assumed that a given crop always has the same treatment 
i.e. that the periods of the crop and treatment cycles are the same. 
Thus we might have an experiment in which treatment cycle A was the 
application of farmyard manure to the potato crop in a potatoes- 
barley-wheat rotation, treatment cycle B being the control (no farmyard 
manure). The barley and wheat crops will then measure the first and 
second year residual effects of the farmyard manure. 

As all the crops are different each phase of the rotation will normally 
be grouped in a separate series of three blocks. Table 2 shows the 6 


TABLE 2. CROP TREATMENT SEQUENCES IN A THREE COURSE ROTATION WITH 
ALL CROPS DIFFERENT 


(Crops: I’, Q, R; treatment eyeles: A; , Az, A3; Bi, Bz, Bs) 


Series: I II III 
Blocks: 1,2,3 4, 5,6 7, 8,9 
Sequence: 1 2 3 4 5, 6 
Preliminary | -1 Q A, Q B, PA, PB, RAs; RB; 
years \ 0 R As R B; Q Ag Q B, F Al P B, 


[experimental 3 R A; RB; Q Az Q Be PA, PB, 
years 


| 
ih 
A | 1 A, P B, R R B; Q A: Q B, 
| 2 | QA: QB | PA, PR | RAs RB; 
an 5 QA, QB, PA, PB, RAs; RBs 
| RAs RBs | QA. Qh | PA PR 
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crop-treatment sequences of such an experiment over the two preliminary 
years and six experimental years. Each of the sequences will be repli- 
cated three times, there being 9 blocks in all of two plots each. The 
yields for the six experimental years for crop P will consist of three 
replicates of Table 3, where a, , 6, , a2 , b2 etc. represent the yields of 


TABLE 3. YIELDS OF ONE REPLICATE OF CROP P 
FROM THE ROTATION EXPERIMENT OF TABLE 2 


Series 
I II Ill 
Year 1 2 3 4 5 6 

1 ay by 

2 a2 

3 a3 bs 
4 as by 

5 as bs 

6 a6 be 


crop P in years 1, 2 etc. and in treatment cycles A and B (suffices 
being used to denote years). The results for the other two crops will be 
similar. 

The analysis in this case is simple. The yields of crop P in series 
I will constitute a 2 X 2 X 3 table of treatments X years (1 and 4) X 
replicates (blocks). The degrees of freedom of the analysis of variance 
of series I can therefore be partitioned as follows: 


Treatments (7') 
Years (Y) 
Blocks (B) 


YXB 
xX B 


Since the treatments are repeated on the same plots the components 
T, B, and T X B which are derived from plot totals over the two years 
constitute the plot total part of the analysis, the remaining terms the 
plot X years part. 


The contrasts between the totals of series I, II and III are estim. .es 
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of year differences (though subject to greater errors than the other 
components of years). The three single degrees of freedom for treat- 
ments from the three partial analyses combine into one for treatments 
and two for treatments X series (i.e. treatments X years). The 
combined analysis of the three series for one crop is therefore as shown 
in Table 4. 


TABLE 4. ANALYSIS OF VARIANCE 


D.F. MS. Expectation 


Series (Years) 2 

Treatments (7') 1 

Plot Treatments X Series = 7’ K Y 2 
Totals Blocks (B) 6 
6 


T X B (plot error) E, 202 + 2 
15 
Years (Y) 3 
Plots X 7x Y 3 
Years |YXB 6 
T X Y X B (plot X year error) 6 Ey o2 
18 
35 


The 7 X B component with 6 d.f. provides an estimate of error 
(which we may term plot error) for comparisons based on plot totals 
over all years, and the JT X Y X B component, also with 6 d.f., provides 
an estimate (which we may term plot X year error) for comparisons not 
involving differences between plots. The estimate of the error variance 
of the difference between A and B averaged over all years, for example, 
will be }£, , and that of the change in this difference between the first 
three years and the second three years will be 4 F, . 

If the errors of the plot yields are regarded as made up of two com- 
ponents, one independent from year to year with variance o2 and the 
other constant over all years with variance o, , the expectation of E, 
is o, and of EL, is 20, + o2. (The coefficient of o; is given by the 
number of yearly yields, two in this case, entering into each plot total.) 

The sum of the sums of squares for the two estimates of error will 
be equal to the sum of the sums of squares for error in the analyses of 
the results of each year separately. There will be 2 d.f. for error in 
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each year giving 12 d_f. in all, and the expectation of each mean square 
will be o2 + o2 , agreeing with the combined expectation of T * Y 
and T X Y X B. 

If results are available for a larger multiple of three years the first 
part of the analysis will be unaltered except for change of divisors, but 
the degrees of freedom in the second part will be increased. With 12 
years, for example, all these degrees of freedom will be multiplied by 
3, and the expectation of /, will be 4 o; +o; . The analysis is also 
easily extended to an experiment containing more than two treatment 
sequences. With ¢ treatment sequences a factor ¢ — 1 will be introduced 
into the degrees of freedom of all components involving T. The altera- 
tions for change of rotation length are equally simple. 

In experiments extending over a long period of time it is sometimes 
advisable, when trends are being considered, to subdivide the plots X 
years part of the analysis further, into say linear trend and remainder, 
so as to take account of possible differential trends of plots within blocks. 
This point is discussed by Cochran (1939). and Patterson (1953) and 
need not detain us here. 

It will be noted that the interaction of years with blocks, Y X B, 
has not been included in the estimate of the plot X year component of 
error. This is as it should be, since the different blocks may well exhibit 
year to year differences: such differences will not enter into the treatment 
comparisons since these are all compounded of differences within blocks 
in the same years. It is the elimination of this component which 
complicates the analysis of the rice-pasture experiment which is the 
subject of our investigation. 


A rotation in which the same crop occurs more than once 


If the same crop occurs more than once in a rotation comparisons 
can only be made with precision between the yields of the crop at the 
different stages of the rotation if the plots carrying them are arranged 
together in blocks. This requires that some at least of the different 
phases of the rotation must be grouped together in the same blocks. 
Whether all phases have to be so grouped depends on the crop rotation. 
In the four course rotation potatoes, wheat, potatoes, barley, for 
example, phases 1 and 3 must be grouped together, as must phases 2 
and 4. In the rotation potatoes, potatoes, bariey, wheat all four phases 
must be grouped together. 

Whether comparisons between the same crop at different stages of 
a rotation are required will depend on the nature of the experiment. 
Such comparisons are usually of some interest, but may not be judged 
of sufficient value to compensate for the loss of accuracy resulting from 
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the larger blocks and the inconvenience and other troubles arising from 
the growing of more than one crop in the same block. We shall not 
discuss the possibilities here, as the problems are essentially those of 
design. 

When the different phases of a rotation are grouped together new 
problems of analysis arise. We will first consider the simple example 
of a three-course rotation with two treatment cycles A and B as before, 
but with one crop repeated in two consecutive years e.g. with the 
rotation of crops P.P.R. instead of the rotation P.Q.R. of Table 1. Each 
replicate of the three phases of the rotation will be taken to be arranged 
together in a block instead of in three separate blocks. The yields of 
one replicate of crop P will now be as shown in Table 5. 


TABLE 5. YIELDS OF THE DUPLICATED CROP IN ONE REPLICATE OF A THREE 
COURSE ROTATION WITH TWO CROPS THE SAME 


Sequence 
= 
Year 1 | 2 3 { ) 6 

| 

—— 
1 bd; at 
2 as bs as by 
3 | af bf a3 | bs 
4 | as bs at by 
5 as bf as bs 
6 | | bf dé bg 

| | 


Here aj , bi , etc. represent the yields of the second P crop under 
treatment cycles A and B. The contrast between 3 (@ + 6) and 
3(a’ + 6’) is a crop-phase contrast. 3(@ — 6 — a’ + 6’) will then repre- 
sent the treatment X crop-phase interaction. 

The comparisons of chief interest will be between the means over 
all years of a, a’, b and b’. The difference between 3(@ + a’) and 
3(b6 + 6’) is derived from plot (i.e. sequence X replicate) totals over 
all years. The differences between d and a’, and between 6 and b’ do 
not involve plot differences. 

If the same crop were grown in all years the interaction component 
of the (6 X 3) plot total X replicate table would represent plot error. 
In the present case, however, the plot totals over the six years will be 
affected by year differences. Sequences 1 and 2 contain years 1, 2, 4 
and 5, sequences 3 and 4 years 2, 3, 5 and 6, etc. Hence the interaction 
component of the plot total X replicate table will contain components 
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of block X year interactions, which as pointed out above, require to 
be eliminated from error. 

A partial analysis is easily effected. If plots 1 and 2, and the corre- 
sponding plots in the other replicates, are taken, we obtain a2 K 4 X 3 
table of treatments X years X replicates, phases being completely 
confounded with years. The interaction of the (2 X 3) plot total X 
replicates table gives 2 d.f. for plot error, and the three factor inter- 
action gives 6 d.f. for plot X year error. Similar analyses can be made 
for plots 3 and 4 and plots 5 and 6. Combining the results we obtain 
an estimate of plot error with 6 d.f. and an estimate of plot X year error 
with 18 

In experiments with a fair number of treatments this partial analysis 
may be judged sufficient to give adequate estimates of error, but it is of 
interest to consider how the remaining degrees of freedom for error can 
be recovered. In experiments containing rotations of varying length 
or with variations in rotations of the same length, without other treat- 
ments, as in the rice-pasture experiment under consideration, recovery 
of these degrees of freedom is essential. 

The number of degrees of freedom for error from a single year’s 
results will be 6. Hence there will be 36 d.f. in all for error. The full 
table of plot totals will consist of a 6 X 3 table of sequences X blocks. 
This indicates that the total number of degrees of freedom for plot 
error is 5 X 2 = 10. Hence there will be 26 d.f. for plot & year error. 
We therefore require an additional 4 d.f. for plot error and an as anes 
8 d.f. for plot year error. 

If c, is written for a, + 6, , c{ for af + b{ etc. the yields of one 
replicate can be written as shown in Table 6. 


TABLE 6. COMBINED YIELDS OF ONE REPLICATE 


Sequences 

Year 1 and 2 3 and 4 5 and 6 Total 

1 C1 cf 

ch C2 Y2 

3 cs C3 

4 C4 cf Y, 

5 cg C5 Ys 


| 

— 

| 

| | 
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Only the differences within columns of the c’s and e’’s, which have 
appeared as year differences, have so far been taken into account. We 
shall now consider the full analysis of this table. 

There are 11 d.f. in all which can be partitioned into 


Columns (plot pair totals) 2 
Years 5 
Columns X Years 


Columns and years are not, however, orthogonal. To overcome this 
we may replace the column totals P,,, etc. by the quantities Q, , Q. , Q; , 
given by 


Q, = 2fe, + 3s +3) +e, +5 + 3(Co + 


etc., where Y,,, represents the total Y, + Y, + Y,+ Y; of the year X 
replicate totals for the years in which sequence 1 (or 2) carries the crop. 
All years are equally represented in these quantities, which are therefore 
orthogonal with years. They also represent differences between multi- 
ples of plot totals. Consequently if the quantities are calculated for all 
three replicates the interaction of the resultant 3 X 3 table of Q’s will 
give 4 d.f. for plot error. 

The divisor for the squares of the Q’s, and the expectation of the 
error in terms of o; and o, , remain to be determined. 

For the divisor we have 


— Q = 2fe. — 3. + 3; +05) + 
= (c, —c) + 2e6 +; + 


Since each c is a total of two plots the divisor of (Q, — Q.)’ is 
2{17 X 8+ 2? X 4} = 48. The divisor of each Q’ is therefore 24. 

To evaluate the expectation in terms of o; and ¢. we replace all 
yields of plot 1 by p, , etc. We then have 


Q, — Q. = 6p, + 6p. — 6p; — bp, 


Hence the expectation of the o; component of (Q, — Q.) is 4 X 6° X 
= 1440, . The expectation of the corresponding error mean square 
is therefore 30, + o.. If years and plots had been orthogonal the 
expectation would have been 403 + o2 . 

Alternatively the sum of the coefficients of the p® terms in Q] + 
+ Q — 4(Q, + Q + Q;)’ can be calculated. 

If these 4 d.f. for plot error and the 6 d.f. already obtained are 
combined by adding the sums of squares the expectation of the re- 
sultant mean square (10 d.f.) will be 
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10 


6 


2 2) _ 3 642 
10 (40, + o,) = 3.60, + 0. 


(30; + on) +- 
This method of combination is only fully efficient if o; is small relative 
to o.. Otherwise greater weight should theoretically be given to the 
mean square with larger expectation. Even in the extreme case when 
o, is very large relative to ¢; , however, the loss of information with the 
above weighting is quite trivial. Nevertheless, as will be seen later, the 
combined estimate is by no means equivalent to an estimate with 10 
d.f. and expectation 403 + o2 . 

The analysis of the three replicates of the c’s can now be completed 
from the table of the totals of the c’s over all replicates. It will take 


the form shown in Table 7, where the phase contrast ¢ — é’ is represented 
by P. 


TABLE 7. ANALYSIS OF VARIANCE OF TABLE 5 


1. Years (ignoring treatments) 5 | From the block 
2. Blocks (replicates) 2>X year totals 
3. Y X B (ignoring treatments) 10 } table 

4. Q’s: totals over blocks (P X Y) 2 ait 

5. Q’s: interaction with blocks (plot error) 4)~ 

6. Remainder: totals over blocks (P, 1; P X Y, 3) 4| By diff 

7. Remainder: interaction with blocks (plot X year error) 8{ Serene 


w 


The sums of squares for the last two items are obtained by differences. 
The total sum of squares for the table of the totals of the c’s ove: all 
replicates equals the total of the sums of squares for items 1, 4 and 6. 
The total of the sums of squares for all seven items is equal to the sum 
of squares for the 36 c’s. Item 7 gives the remaining 8 d.f. for plot X 
year error. 

It will be seen that the key to this analysis is given by the expressions 
Q which make the plot pair (column) totals orthogonal with years. In 
this example these expressions can be written down by inspection, but 
where this cannot be done they can be obtained by fitting constants 
for years (rows) and columns by the method of least squares. In the 
next section the method will be illustrated by application to the rice- 
pasture experiment. 

The full analysis can now be completed. If required T, P, T X P 
and their interactions with years can be separately exhibited. It then 
takes the form shown in Table 8. 
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TABLE 8. FULL ANALYSIS OF VARIANCE 


IDF. Expectation D.F.| Expectation 
T 1 | 402 + o2 | Years (ignoring treatments) | 5 
P 1 | Blocks 2 
1 | o2 | Y X B (ignoring treatments)! 10 
| 5 | 0402 + | Plot error 10 | 3.665 + of 
5 | 12 +o, | Plot X years error | 26 
TXPXY}.5 71 


T, P,T X P,T X ¥ and T X P X Y are computed in the ordinary 
manner. P X Y is obtained by subtraction of P from items 4 and 6 of 
Table 6. The errors of the various comparisons and the expectations of 
mean squares such as T X Y can be evaluated in terms of o; and o; 
by replacing yields by plot constants in the manner already indicated. 
The mean expectation of the 18 d.f. for treatments and treatments X 
years checks to ¢, + o; , as it should. 

It should be noted the plot error mean square no longer gives the 
error of differences based entirely on plot totals. Thus the error 
variance of the treatment difference }(a@ + a’) — 3(b + is + 02). 
The estimate of this will be 2(,°h, — }£,). The accuracy of this 
estimate of error, relative to the direct estimate based on the 6 d.f. of 
the partial analysis is about equivalent to an extra 2 d.f. when oa; is 
small relative to o; , but equivalent to nearly an extra 4 d.f. when a; is 
large relative to ox . 

The analysis of the other crop of the rotation will be similar to the 
analysis of Table 4. There will now only be 2 d.f. for blocks, the re- 
maining 4 d.f. representing year X block interactions. There will still 
only be 6 d.f. for plot X year error, since the other 4 d.f. will be com- 
pletely confounded with year X block interactions. 


Orthogonal partition of the degrees of freedom in the rice-pasture expeiment. 


In the rice-pasture experiment orthogonal Q functions for vhe 11 
crop sequences cannot be written down by inspection. Fitting of 
constants must therefore be resorted to. The required functions can 
be obtained from the equations appropriate to a single replicate. 

The yields and year and plot totals of one replicate are shown in 
Table 1. Let the year constants be y, , yz --+ yi2, and the plot constants 
P:,P2,°** Pu. The normal equations are then: 
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Denote the sum of the Y’s for the years in which plot 1 carries rice 


4y, + pi + ps + Pr + Ds 
Ay, + po + ps + Ds + Do 


Y, 
Y, 


P, 


P, 


Ps 


by Yay so that Ya = + + + etc. 


+ pu = Ps-1. , With similar totals for the P’s. 


Eliminating the y’s in the p equation by means of the y equations 


we then have: 


12p, 
12p, 
12ps 
18p, 
18ps 
18p5 
IDs 
Ips 


— 2Ps,7 — Ds-11 
2Ps.6 — 2Ps,7 — Pa-i1 
2Ps.6 — — 
— 3Ps,7 — 3ps 
2pi-s — — 3Ds 
Pi-s — 3p, — 3p; 
Pi-s — 3p, — 3ps 
Pi-s — 3ps — 3p 
Pi-s — — 3p; 


3Pr0 


= 


— 3pu 


= 4P, 
= 
= 4P, 


= 4Py, — 


= 


Yay 
Y (2) 
Yu 
Ys) 
Yw) 
Yin 
Ys) 
Yio) 
Y ao 


} (il) 


We may now take the following sums of these equations: 
(1) + 2) + (3) 

12p,.3 — — Gps-7 — = 

3(4) + 3(6) + (5) + (7) 

— — = + 4Ps.2 — 


Further put 
Pi + Po + Ps = Di-a, Pi t+ Do = Dao, Ps + Pz = Ds.7 » Ps + Po + Pro 
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(1) 
(2) 
(3) 
(4) 
(5) 
(6) 
(7) 
(8) 
(9) 
(10) 
(11) 


4p, +yt+ + Y7 + = : 
+ yi + Ys + Yo 
=> 4P 1 
=4P, 
= 4P, — 
- = 4P, 
= 4P, — 
| | 
4R 
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3(5) + 3(7) + (4) + © 
36p5.7 — 16p,.3 — 12ps_,, = 12P;.7 + 4Py., — 4R 
(8) + (9) + (10) + (11) 
9pe11 — 41-3 — 3ps.6 — 3ps.7 = 4Ps1, —R 
It is easily verified that these equations are satisfied by 


(12) 
= 
Pe-11 = $P 3-11) 


This is a consequence of the fact that the sums of plots 1-3, 4 and 6, 
5 and 7, and 8-11 are all orthogonal with years. 

The values of the separate p’s are then obtained by substitution in 
equations (1)—(11), and elimination of ps — p,, from equations (4)—(7) 
by means of equations (8)—(11). This gives 


36p, = 17P, — Pro + + — 5Y + + Pi-s 
24 4,6 2* 5.7 


The total sum of squares accounted for by fitting the constants 
might be calculated from the expression S(yY) + S(pP), the values 
of the year constants y being determined by substitution of the p values 
in the original normal equations, or alternatively from the expression 
+ S(pP’), when = 4P, — ete. It is more satisfactory, 
however, to construct a set of orthogonal Q functions which partitions 
the sum of squares directly. 

Such a set of functions can be constructed by taking appropriate 
groups of constants in turn in any desired order, and obtaining ex- 
pressions for these in terms of the Y’s and P’s after eliminating all 
groups of constants already fitted. If certain conditions of symmetry 
are satisfied the sum of squares attributable to a group of constants will 
be calculable from the sum of the squares of the deviations of the values 
of the corresponding Q functions in the group. 
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As a first step we may replace the p’s by fi-3 , Ps.5 , Ps.7 and Ps-11 , 
5p, , , Spy, where = 3Pi-s, Pais = ete. and dp, = 
Di — = Ps — , etc. We may then take the following 
groups of constants in the order shown: 

Yi, Y2,°°° 

Pi-s » y Ps-11 
Opi , , Ops 

. Ops , , 
ops and Ops op; 


no 


The ordinary expressions of the analysis of variance give the sums 
of squares accounted for by groups 1 and 2, which are orthogonal. 

The values of 6p, , dp. , dp; before fitting dp, to 6p,, are given by 
replacing ps, by j4,. etc. in equations (1)—-(3). It will be seen that the 
equations (13) already obtained are unaltered, and the corresponding 
sum of squares will be given by 


where Q, = 4P, — Y,, etc., dev’ indicates the sum of the squares of 
the deviations, and d is a divisor to be determined. 
We have, by ordinary least squares procedure, 


V(é6t, — = (Cur — + 


where c,, and ¢,, are the values of 6t, and 6¢, when P, = + 3, P, = 
P, = Pas = = Ps-1, = 0 and all Y = 0, and (= 
and c,. are the values of é¢, and 6¢, when P, = 2, P, = P; = —}, P44 = 
P;,; = Ps_., = 0, and all Y = 0. From equations (13) 


+ 1 
12 dp, = 4P, — Yo —3Pistgk 


Thus 


= 
| 
aL 
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Hence 


« 


~ 


2 
o 


— 6) = V (Q, — = 


Hence 


The values for dps to ép,, before fitting dp, to dp, are given by equa- 
tions (8)-(11), replacing p, and p, by j,,, and p, and p, by p;.; . Hence 


Ops 4P, Y¢s) + iP,.s + + 


By putting Ps = +2, Py = Py = P,, = —} and all other P and all 
Y = 0 we obtain css = +4, Coo = —4, and therefore css — 2cgo + 
Coo = §. Hence the divisor is 36. 


Finally the values for 5p, to dp, after fitting all other constants are 
given by equations (14). We have 


12( dp, 5ps) 3(P, — PF) 
+ Po t+ Py — Pro — Pu Ys) Yo) + Yao + Yaw 
By putting P, = +3, P; = —} and all other P and all Y = 0 we obtain 


Cag — Cog = 
and therefore cy, — 2css + Cos = 3. Hence the divisor of 144 (dp, — 
is 72. 

We may therefore summarize the analysis of a single replicate as 
follows: 

In addition to the year totals, Y, — Y,,. , and the totals of the four 
groups of plots orthogonal with years , , , calculate 
the following 11 quantities Q, — Q,, . 


For plots 1-3: Q, = 4P, — Yq), ete. Total: 4P,_, —R 

47: Q, = + Ps + Py — — Yo) ete. 
Total: 3P,.6 + 3Ps,7 + 2Ps_1. — 2R 

» $-11: Qs = 4P, — , etc. Total: 4P;_,;, — R 


The expressions for the sums of squares are shown in Table 9. 

The reason why the groups were taken in the chosen order will now 
be clear. If 6p, — 6p,, had been taken after 5p, — 5p, equations (15) 
would have been used for determining cys , Cs9 , etc. In this case cs, 
would not be equal to cs 49 , ete., i.e. not all the covariances between 
éps — 6p,, would be equal. Consequently the expression for the sum 
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TABLE 9. ANALYSIS OF VARIANCE: EXPRESSIONS FOR SUMS OF SQUARES 


D.F. S.S. 


1 
. . 1 
Plot Within groups 1-3 2 18 dev'(Q, , , Qs) 
Totals) groups 4,6 and 5,7 2 — + (Q; — Q:)°} 
group 8-11 3 dev’ (Qs Qs Qio 
| Total 10 
243 
Remainder (plots X years) = 
of squares for 6p, — 6p,, would be more complicated, without any 


compensating simplification in the sum of squares for 6p, — 6p; . 

The divisors in the above sums of squares may be checked by setting 
out the coefficients of the differences of pairs of constants and summing 
their squares. Thus 9(é6p, — 6p) has the coefficients shown in Table 10, 
with sum of squares 72. 


TABLE 10. COEFFICIENTS OF 9(éps — 6ps) 


Year 1 2 3 4 5 6 7 8 9 10 11 Total 

Total o © OO 0 +8 0 -3 +9 -9 0 0 0 


Similarly the orthogonality of the different components can be 
checked by verifying that the sum of the products of the corresponding 
coefficients of the differences of each pair of constants is zero. 

An alternative numerical check of the divisors and of orthogonality 
can be made by assuming values of the year and plot constants, building 
up the plot yields from them, and carrying out the analysis. If the 


Fotis 
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expressions are correct the residual sum of squares (26 d.f.) should be 
zero. This provides a useful overall check and a test of the numerical 
procedure, though it should be remembered that a single numerical test 
does not necessarily detect all errors. Taking y, = 1, yz = 2, --+ Yio = 
12, and p, = 1, p, = 2, p,, = 11, the yields of plot 1 are 2, 5, 8, 11, 
those of plot 2 are 4, 7, 10, 13, ete. and the analysis shown in Table 11 
is obtained. 
TABLE 11. NUMERICAL CHECK 


D.F. S.S 
Years 11 659.75 
Orthogonal groups 3 344.25 
Qi, Qe, Qs 2 6 
Qs — Qs — 2 16 
Qs, Qs, Qo, Qu 3 7.25 
21 1033.25 


The total agrees with the total sum of squares of deviations of all yields 
(47 d.f.) 
Scheme of analysis of the rice-pasture experiment 


We may now complete the scheme of analysis for the rice-pasture 
experiment. Table 12 shows the general scheme, which has the same 
pattern as Table 8. 


TABLE 12, ANALYSIS OF VARIANCE OF THE RICE-PASTURE EXPERIMENT 


D.F. 

( Crop sequences 10 

Plot Blocks 2 

Totals Plot error 20 
| 

| Total 32 

Years 11 

Crop sequences X years 26 

Plots X | Blocks X years 22 

Years Plot X year error 52 


Total 
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Tables corresponding to Table 1 must be prepared for each replicate. 
and for the total of all three replicates. The sums of squares for blocks, 
years and blocks X years are computed in the ordinary manner from 
the marginal totals for years. For crop sequences and plot error the 
quantities P,_, , Ps.6, Ps.2, and Q, — Qi, must be calculated for 
each replicate and for the total of all replicates. These are set out in 
four tables, namely that for the P totals and those for Q, to Q; , Q, toQ; , 
and Q, to Q,, . Columns of the differences Q, — Q, and Q; — Q; are 
also required. The marginal totals for all replicates (with an additional 
factor 3 in the divisors) give the components of the sum of squares for 
crop sequences (10 d.f.) while the interaction components (which for 
Q, — Q,and Q; — Q; are merely the sums of the squares of the deviations) 
give the components of the sum of squares for plot error (20 d_f.). 

The sum of squares for crop sequences X years (26 d.f.) is then 
obtained by subtracting the sums of squares for years and for crop 
sequences from the total sum of squares (47 d.f.) for the crop sequence 
X years table (totals over all replicates). Likewise the sum of squares 
for plot X year error (52 d.f.) is obtained by subtracting all other items 
from the sum of squares (143 d.f.) for the whole experiment. 

It remains to determine the expectations of the various components 
of the plot error sums of squares in terms of o, and o; . The procedure 
already explained of evaluating the sums of squares of the coefficients 
of the plot constants is followed. The group 8-11 component can be 
broken down into 3 orthogonal components. We may take the ortho- 
gonal components 9(éps — 4po), 9(6pio — and 9(éps + Spy — 


TABLE 13. EXPECTATIONS OF PLOT-ERROR MEAN SQUARES 


D.F. Expectation 
Pia Pez 4 60,+ 


Orthogonal groups) P;,; — 2 
+ Ps; — — Psu 2 

Within group 1-3 4 36+ 0, 
groups 4, 6 and 5,7 4 

group 8-11 [dps + — — Spi 2 Bota. 
Total 20 3.60,+ 0: 
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— Theo, component of 81(ép. — 6p)”, for example, obtained 
from the total line of Table 8, is (3° + 3° + 9° + 9°) o = 180.0%. Since 
the divisor is 72 the expectation of the corresponding sum of squares 
(1 df.) is 36, + o; . A similar procedure can be followed for the 
orthogonal groups component. The complete set of expectations is 
shown in Table 13. 

The expectation for the whole 72 d.f. for plot and plot X year error 
is therefore 


| D.F. | M.S. Expectation 
Plot error 20 | E, 3.6 o. +2 
Plot X year error | 82 | # o2 
Total 72 + o2 


This is as it should be, since the 72 d.f. is made up of 6 d.f. for error for 
each year separately, each with expectation o, + o3. The expectation 
for plot error could of course be obtained by utilising this fact, but 
direct evaluation provides a useful check and exhibits the structure of 
the analysis. 

In the above analysis the degrees of freedom for crop sequences and 
for crop sequences X years have been partitioned in a manner suitable 
for the separation of plot error and plot X year error. They can also 
be partitioned so as to isolate the various contrasts which are of ex- 
perimental interest. The main contrasts are those between rotations 
A and C and the two phases B and R’ of B. The contrasts between A 
and C and the mean of B and B’ are part of the 10 d.f. for crop sequences. 
The sum of squares (2 d.f.) is given by 


144 
where 7,_; represents the total over the three replicates of P_; , etc., 
and T the grand total. The corresponding mean square is not, however, 
directly comparable with the plot error mean square, since its expecta- 
tion, from the results already given, is }(} + 42) 0, +o. = 33 
It must therefore be compared with $3 EF, — 7, E,, . 

The contrast between B and B’ is wholly within plots and is therefore 
subject to plot X year error, the error variance being ;', I’, , and the 
sum of squares (1 d.f. from crop sequences and years) ~,;{S(B) — 


S(B’)}’. 
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The remaining 33 d.f. from crop sequences and crop sequences X 
years represent the interactions of rotations A and C' and the two 
phases of rotation B with years. The expectation of the corresponding 
mean square is z's (36 — 33) of + 0. = 340, +o, and it is therefore 
compared with 433 + 

It should also be noted that blocks X years contains components 
of plot differences. 


Analysis of incomplete data 


In long term agricultural experiments an interim analysis of the 
results is often required. The interim data usually lack the balance 
that is attained when the experiment has run its full term. In the type 
of experiment we have just been discussing an exact analysis may then 
prove excessively laborious since the relevant normal equations will no 
longer be readily soluble. On the other hand it is unreasonable to tell 
the experimenter that he must await the completion of the experiment 
before attempting to draw any conclusions from the results. The 
statistician must therefore be prepared to determine how far approxi- 
mate methods of analysis will enable interim reports to be made on 
long-term experiments in an expeditious manner and without undue 
labour. 

In the case of the rice-pasture experiment there do not appear to be 
any approximate methods which are very satisfactory. As an example 
of the problems involved we may consider the analysis of 6 years’ data. 
We may adopt the same notation as previously except that P, etc. and 
T, etc. now represent totals over 6 years instead of 12. In the exact 
analysis on the lines already laid down, the degrees of freedom would 
partition in the manner shown in Table 14, which corresponds to 
Table 9. 


TABLE 14. ANALYSIS OF VARIANCE OF 6 YEARS’ DATA 


D.F. 

| Crop sequences 10 

Plot Blocks 2 
Totals Plot error 20 
Years 5 

(Crop sequences X years 8 

Plots X Blocks X years 10 
years { Plot X year error 16 


| 
| 
71 
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An estimate of plot error could be built up from the orthogonal 
groups and the Q functions already obtained. The components derived 
from the Q functions are now, however, no longer fully orthogonal 
amongst themselves (though they are still, as is essential, orthogonal 
with years). This will not seriously affect the estimate of plot error— 
it has an analogous effect to the inclusion of some degrees of freedom 
more than once in an estimate of error—but the sum of squares so 
obtained cannot now be subtracted from the total sum of squares for 
error to give the plot ” year error without serious risk of disturbance 
due to non-orthogonality. 

The simplest procedure, therefore, is to throw the block X year 
interaction into the estimates of error. The interaction of the table of 
treatment sequence totals X blocks (appropriate allowance being made 
for differing numbers of yields in the totals) will then give an estimate 
of plot error (containing components of block X year interaction) with 
20 d.f. The total of the interaction components of the 11 tables of 
treatment sequences X replicates (the last two of which contribute 
nothing) will give an estimate of plot X year error (also containing 
block X year interaction) with 26 d.f. The expectation in terms of 
g, and o. of the plot error can be evaluated by the procedure already 
adopted. 

This method of estimation of error will give estimates which are 
too large if block X year interactions are substantial. This point can 
be examined by comparing the mean square for block X year inter- 
actions, calculated from the table of block X year totals, with its error 
expectation in terms of o; and oc. . If it is considered that the inter- 
action should be eliminated then a full least square analysis, with 
inversion of the matrix, will have to be undertaken. Construction of 
orthogonal functions, which will necessarily be complicated, will not 
be worth while. 
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THE DERIVATION OF JOINT DISTRIBUTION AND 
CORRELATION BETWEEN RELATIVES BY THE 
USE OF STOCHASTIC MATRICES 


C.C. Li anp Lovis Sacks 


Department of Biostatistics, Graduate School of Public Health, and Department of 
Mathematics, University of Pittsburgh, Pittsburgh, Pa. 


The absolute frequencies of the various genotypic parent-offspring 
and sib-pair combinations with respect to one pair of genes, autosomal 
or sex-linked, are well known to geneticists, for they are of practical 
value in certain types of studies in human heredity. These frequencies 
are, however, obtained by a rather long process, even for the simple 
case of full sibs. The procedure of obtaining the frequencies of uncle- 
nephew or first cousin combinations is entirely too tedious even with 
the help of matrix notations (Hogben, 1933). The purpose of the 
present communication is three-fold. The first is to give a simple 
procedure of finding the frequencies of the various genotype combi- 
nations of near relatives by using matrices of conditional probabilities. 
The second purpose is to express such matrices of conditional probabili- 
ties of relatives in the form of a linear function of some basic matrices. 
The third is to deduce the correlations between relatives from such 
linear functions of the basic matrices. The meaning of these statements 
will be made clear in later sections. 

For many years there existed two apparently very different methods 
of obtaining the genotypic correlations between relatives. One is a 
straightforward but long procedure by which the frequencies of the 
various combinations of the stated relatives in the general population 
are first found and then the correlation is calculated from such a “‘corre- 
lation table’. On the other hand, they may be obtained by the method 
of path coefficients, developed by Wright (1921 and later). The latter 
method gives the required correlation coefficient almost instantly once 
the relationship is specified, but does not give us any information about 
the frequencies of the various combinations of the relatives in the 
population. Moreover, one has to be familiar with the mathematical 
theorems concerning the path coefficients before he can use them to 
derive the required correlations. The method of stochastic process 
and its final reduction to some basic matrices, as described in the 
following pages, will give us both the frequencies of relative-pairs and 
their correlation in a very simple manner. It is hoped that the present 
method will bridge the gap between the two existing procedures. 
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AUTOSOMAL GENES 


1. Instead of directly deriving the absolute joint distribution with 
regard to the genotypes of certain relative pairs (e.g. parent-offspring, 
uncle-nephew, first-cousins, etc.), we shall deal with the conditional 
probabilities that a relative should be of a certain genotype when the 
other relative’s genotype is given. Thus, when an uncle is given to be 
AA, we are to find the probabilities that his nephew should be AA, 
Aa, aa, repectively. These three probabilities will of course add up to 
unity. Similar probabilities may be found for Aa and aa uncles. Thus, 
we obtain the following array of probabilities (¢): 


Nephew 
AA Aa aa 
Condition: Uncle AA bes 
Condition: Uncle Aa bse 
Condition: Uncle aa hes i. too 


where the sum of the three probabilities in each row is unity. Such an 
array of probabilities is known as the matrix of “transition probabilities”’, 
from uncle to nephew in-this case. Once such a matrix is obtained, the 
absolute frequencies of the various genotypic combinations of uncle and 
nephew in the general random mating population may be obtained 
immediately by multiplying the uncle conditions by their respective 
“initial” probabilities; that is, multiplying the first row by p’, the 
second row by 2pq and the third row by q’, where p is the frequency of 
gene A in the population (¢ = 1 — p). Therefore, we shall largely deal 
with the transition matrices of relatives in the following. 

2. In order to facilitate the derivation of transition matrices for 
various types of relatives it is convenient to introduce another con- 
sideration besides the genotypes. Let X be a gene of the locus under 
consideration. It may be A ora. We say that it exists in two states. 
If we select a gene from each of two unrelated individuals and they 
happen to be both A or both a, we say that these two genes are alike 
in state. On the other hand, if a parent has an X-gene which is trans- 
mitted to his child, we say that these two X-genes (one drived from the 
other) are identical by descent. Genes of the latter kind are necessarily 
alike in state, barring mutation, but the reverse is not true. The 
distinction between the two ways in which two alleles may be alike 
has been pointed out by a number of geneticists (for example, Cotterman 
1941; Malécot 1948; Crow 1954). To trace the origin of a gene, we may 
label each independent gene by a subscript (a serial number for identi- 
fication). Thus, in a random mating population, the two parents may be 
represented by X,X, X X;X, , whatever their genotypes. 
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3. A pair of relatives may have both, one, or none genes identical 
by descent, depending upon the type of their relationship. The matrices 
of transition probabilities (with regard to genotypes) for the three 
cases are respectively: 


p p pq ¢ 
1 O|, T= |4p 4 O= Bg QI. 
001 0 pq p pq ¢ 


The truth of J and 0 is obvious. If two relatives share both genes 
identical by descent, they are necessarily of the same genotype. On 
the other hand, if they share no gene identical by descent, the probabili- 
ties that the other relative’s being AA, Aa, aa are p’, 2pq, g’, whatever 
the genotype of the given relative. The truth of 7 is also easily seen. 
For instance, if one relative is given to be AA and the other relative 
shares one gene identical by descent, then the probability that the 
latter should be Aa is q, which is the probability that the other (non- 
identical) gene is a. 

4. With the above three basic transition matrices, it is easy to 
derive the conditional probabilities for any specified relatives. Let us 
consider first some of the simple wnilineal relatives. Now, a pair of 
parent and offspring necessarily share one gene identical by descent. 
Therefore, the transition probabilities from parent to offspring are 
simply the elements of 7. This may be readily verified by the long 
classical method of considering all the possible mating types in the 
population. It should be noted that T also gives the transition probabili- 
ties from offspring to parent. 

Next, consider the transition probabilities from grandparent to 
grandchild. If the grandparent is given to be AA, the probabilities of 
the genotypes of the parent are the elements of the first row of T. For 
each genotype of the parent, the probabilities that the child should be 
AA are the elements of the first column of 7. Thus, the total con- 
ditional probability that the grandchild should be AA is the sum of the 


products of the corresponding elements of the first row and of the first 
column of 7’; viz. 


Pp 
(p,q, 0)|3p| = p + = 4p + 
0 


Similar considerations show that the conditional probabilities for the 
various genotypes of grandchild with one grandparent specified are 
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given by the elements of the product matrix 7-T = T’. Labeling and 
tracing the origin of genes shows that a grandparent-grandchild pair 
has an equal chance to share one gene identical by descent and none at 
all. Therefore, 


T° = 3T + 30 


as may be readily verified by direct multiplication of 7 by 7. Since 
T represents the transition matrix from child to parent as well, T” also 
gives the conditional probabilities for half-sibs, who have one parent in 
common. The above relation shows that the tedious multiplication of 
matrices is reduced to the simple operation of addition. Further in- 
vestigation shows that, 


90, 
and in general, 
= (4)"T + (1 — (3)")0 


where n + 1 is the number of generations between the two relatives. 

5. Now, we may turn to bilineal relatives, the most important type 
of which is full sibs. Let the given sib be X¥,X;. Since X, comes from 
one parent and X, comes from the other, his parents may be designated 
by X,X, X X3;X, , in whatever states the genes may exist. Then the 
probability that the other sib should also receive X, from the first 
parent and X; from the second parent is }. That is to say, the prob- 
ability that the other sib should have both genes identical by descent 
is 3. This probability is independent of the genotypes of the sibs, 
independent of the genotypes of their parents and independent of the 
gene frequencies of the general population. Therefore, the three kinds 
of sib-pairs, AA—AA, Aa-—Aa, aa-aa, should have a component of } in 
all populations. In fact, such sib-pairs are equivalent to ‘identical 
twins’ as far as the A-a locus is concerned. 

On the other hand, the probability that the other sib should be 
X,X, , sharing no gene identical by descent with the given X,X; sib, 
is also}. Finally, the probability that the sibs share one gene (X, or X;) 
identical by descent is 4. Hence the transition probabilities for full 
sibs are the elements of the matrix. 


S = 31+ 47 + 40. 


This, again, may be verified by the long classical method of considering 
all possible matings in the general population. 

Another relationship in which the relatives can share both genes 
identical by descent is that of double first cousins, whose parents are 
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members of two sibships. There are six types of sibships in the popula- 
tion, corresponding to the six types of matings. Double cousinship 
can occur in twenty-one distinct ways, since any pair of the six types of 
sibship may be taken. The absolute frequencies of the various double 
cousin pairs in the population may be calculated by random mating such 
pairs of sibships. ‘This is the procedure employed by Fisher (1918) and 
it involves a large amount of algebra. The conditional frequencies of 
double cousins, however, may be obtained in the following simple 
manner. 

We label and trace the eight independent genes of the four grand- 
parents involved in the two families, as illustrated in Fig. 1. When one 
(any one) cousin is given, the probabilities that the other cousin should 
have both, one, or none genes identical by descent are 7, 3%, 15, 
respectively. Hence the matrix of conditional probabilities for double 


1 92 3,4 5,6 1,8 


O O O O 


Given: 1,5 1,5; 1,6; 1,73 1,83 
2,5; 2,6; 2,73 2,8; 
3,5; 3,6; 3,7; 3,85 
4,5; 4,6; 4,75 4,8. 


FIG. 1.) RELATIONSHIP BETWEEN DOUBLE FIRST COUSINS. THE NUMERALS 
ee 8 ARE LABELS OF THE EIGHT GENES OF THE FOUR GRANDPARENTS INVOLVED 
IN THE TWO FAMILIES. THE FOUR PAIRS OF LABELS BELOW THE PAREN7S INDI- 
CATE THEIR POSSIBLE GENETIC CONSTITUTION ACCORDING TO ORIGIN OF GENES, 
ETC. FOR EACH GIVEN COUSIN, THERE ARE 16 POSSIBILITIES FOR THE OTHER 
COUSIN, OF WHICH ONE IS IDENTICAL (1, 5), THREE WITH GENE: BUT WITHOUT 
GENEs , THREE WITH GENE; BUT WITHOUT GENE:—A TOTAL OF SIX POSSIBILITIES 
SHARING ONE IDENTICAL GENE. THE REMAINING NINE GENOTYPES ARE INDE- 
PENDENT OF THE GIVEN CONDITION. 


1,3; 1, ks 5,73 5,8; 
2,3; 2,4; 6,7; 6,8. 
{ 


352 BIOMETRICS, SEPTEMBER 1954 


first cousins is, 
= ig + + 16 


This matrix, .it is to be noted, happens to be the ‘‘square’’ of the S 
matrix for siblings. 

6. Now, we shall consider the unilineal relatives for whom both S 
and 7 matrices are involved. The term “‘uncle-nephew” is used here 
in its broad sense, including uncle-niece, aunt-nephew and aunt-neice 
relations as we are dealing with autosomal genes. When an uncle is 
specified, the conditional probabilities for his nephew (his full sib’s 
child) will be given by the product matrix ST. Conversely, when a 
nephew is given, the conditional probabilities for his uncle (his parent’s 
full sib) will be given by the product matrix 7S. It may be easily 
verified that ST = TS, showing that the uncle-nephew transition 
matrix is also the nephew-uncle matrix. But, the most remarkable 
property is that 


ST =TS =T 


indicating that uncle-nephew relations are the same as those for grand- 
parent-grandchild or half sibs, as established by other methods in 
human genetics. 

Similarly, one’s first cousin is his parent’s full sib’s child; hence the 
conditional probabilities for first cousins are the elements of the matrix 


TST = T° 


which also gives the conditional probabilities for the great-grandchild 
when one great-grandparent is given. In a similar manner other 
relatives may be expressed by a product matrix involving T and S. 

In the light of the results derived above, we see that Feller’s con- 
clusions (1950, pp. 102-103) that ‘“‘brothers are as close to each other 
as grandfather and grandson” (equivalent to saying S = 7”) and that 
the conditional probabilities for a great-grandson or a nephew when a 
great-grandfather or an uncle is given are the same (equivalent to 
saying ST = T*) are both incorrect. These errors are due to his failure 
to use the sibling transition matrix S in his calculations. 


SEX-LINKED GENES 


7. The above method of conditional probabilities may be applied 
to sex-linked characters with only slight modifications. Without loss 
of generality, we let the heterogametic sex be males (A or a) and the 
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homogametic sex be females (AA or Aa or aa) as in the case of mankind. 
There will be four parent-offspring transition matrices (7) according 
as the parent and offspring are males or females. To distinguish them, 
we shall use two subscripts, the first to denote the number of rows and 
the second the number of columns in the matrix of transition probabili- 
ties. Thus, , , . are the matrices for mother-daughter, 
mother-son, father-daughter, father-son, respectively. The matrix 
T3 is the same as 7 for autosomal genes. The remaining three are as 


follows: 
0 
al. q | m= [P | 
01 P 4q 


As before, 0 denotes an independence matrix, with appropriate sub- 
scripts. Thus, 


T's2 = 


2 
032 = 0.3; = i Pa 
pq Pp Pd 


In this notation, the father-son transition matrix T,, = 0.2, indicating 
that the genotypes of father and son are independent. Similar to the 
T’s, there are four sibling transition matrices, according as the given 
sib is a male or female. The following are for sister-sister, sister-brother, 
brother-sister, brother-brother, respectively (subscript 3 = sister, 2= 
brother): 


S33 = 3133 + 47 33 32 + 3032 
So3 = 37 23 + 3023 . Soo = 3102 + 2022 : 


The coefficients in the linear expression of each S give us the nature of 
the siblings. For example, S.. = 3/.. + 30.. means that the two 
brothers either have the same identical gene or have two independent 
genes. 

8. The matrix of conditional probabilities for other relatives may be 
obtained by multiplying the appropriate matrices, just as in the case of 
autosomal genes. But the kinds of relatives for sex-linked genes are 
many times more than the corresponding relatives for autosomal genes 
as each individual involved has to be distinguished in regard to sex. 
It is unnecessary to enumerate all the near relatives. We shall illustrate 
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the method by considering grandparent-grandchild and relatives of the 


“uncle-nephew”’ type. 


each type. 


There are eight different kinds of relations for 


the elder relative is being given. 
Note from the upper part of Table 1 that whenever two males are 


The results are given in Table 1 in which we assume that 


TABLE 1. MATRIX OF CONDITIONAL PROBABILITIES FOR SEX-LINKED GENES 
When the grandparent is given 
grandmother- grandfather- grandmother- grandfather- 
granddaughter granddaughter grandson grandson 
maternal | | T33T2=}T2+ 302 | =}I2+ 
paternal = = O23 T2Tx = Ox = On 
When the Aunt or Uncle is given 
aunt-niece uncle-niece aunt-nephew uncle-nephew 
paternal | | = On = On 


involved in successive generations (i.e. a father-son step) the result is 
an 0 matrix, indicating independence between the two specified relatives. 
But, maternal grandfather-grandson are genotypically related, having 
a probability of 3 of sharing an identical gene. 

Another interesting point to be noted is that: 


1 0 p 
0 
T32T23 = 3 4 = |3p 3¢| =T 


That is, the relationship between paternal grandmother and grand- 
daughter is the same as that between mother and daughter as though 
the intermediate father were non-existent. ‘That, of course, should be 


the case because the father always transmits the same gene he received 
from the grandmother. This product of two transition matrices is also 
equivalent to a “factorization” of the 7 matrix for autosmal genes into , 
two steps: the first step is to specify the probabilities that a gene is 
transmitted from parent to offspring, and the second, the probabilities 
of forming the three genotypes with that one gene given. 
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CORRELATION, GENOTYPIC \ND PHENOPY PIC 


9. The classical method of caleulating the genotypic correlation 
between relatives is first to derive the joint distribution (correlation 
table) of the relatives in the general population. This can be obtained 
by multiplying the rows of our final transition matrix for the relatives 
by the appropriate “initial” frequencies of the given conditions. The 
correlation table thus obtained need not be given here. Suffice to say 
that the results are in complete agreement with those given by Fisher 
(1918) for autosomal genes and those by Hogben (1932) for sex-linked 
genes. 

10. The correlation between relatives, however, need not be calcu- 
lated from their joint distribution in the population. Before we show 
how it may be obtained directly from the components of the transition 
matrix of the stated relatives, we will first make an observation con- 
cerning the additive property of the correlation coefficient under some 
restrictive conditions. It may be shown that if there are N, pairs of 
values in one population with certain mean values, variances, and 
correlation coefficient r, ; and if there are N, pairs of values in another 
population with the same mean values and variances but correlation 
coefficient r, , then the correlation between the NV, + N, pairs of values 
in the pooled population is 


N, N, 
That is, the correlation in the pooled population is the weighted mean 
of the separate correlations. Or, in other words, the pooled correlation 
may be subdivided into two components, corresponding to the con- 
tributions of the component populations. This analysis may obviously 
be extended to any number of component populations as long as they 
have the same mean and variance. 

11. Now, consider the transition matrix for a specified type of 
relatives. Its general expression is, for autosomal genes, 


+¢,T 


where ¢; , Cr , Co are the probabilities that the specified relatives should 
have both, one, or none genes identical by descent (ec; + er + cy) = 1). 
The J-component contributes perfect (unity) correlation between the 
genotypes of the relatives while the 0-component contributes nothing 
to their correlation. If r7 is the correlation for the 7-component (to 
be calculated after converting 7’ into absolute frequencies), then the 


A 
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total correlation between the relatives in the general population will 
simply be 


r=C, + r 


on account of the additive property of correlation coefficient described 
above. Therefore, to calculate the total r, we need only to calculate 
the one basic correlation ry once and for all. Multiplying the three 
rows of T by p’, 2pa , q’, it will be found that, as is well known, that 
rr = 3. Hence, for example, the correlations between 


grandparent-grandchild,r = 0 + (4) = 


2 \2 
double first cousins, 16" i6+ 16 (3 


It should be noted that although these two types of relatives have the 
same correlation coefficient, the joint frequency distribution (corre- 
lation table) of these two types of relatives in the general population are 
quite different; for, in the former case, the transition matrix is }7 + 30; 
in the latter, #7 + 37 + #0. This situation is exactly analogous 
to that of parent-offspring vs. full sibs. 

12. One great advantage of subdividing the correlation between 
relatives into two components, as given in the preceding section, is that 
when A is dominant to a, the above formula for r still holds except that 
rr = 3 is to be replaced by r} = g/(1 + q), the latter being the correla- 
tion in the 7-component population with dominance (Li, 1954, Chap. 
2). For example, the phenotypic correlation between full sibs is 


in agreement with the more conventional expression (1 + 3q)/4(1 + q). 
It assumes the value 5/12 when g = 3, as obtained by Yule (1906, ef. 
Fisher, 1918). The phenotypic correlation between other types of 
relatives may be obtained in a similar manner. 

13. For sex-linked genes the procedure of deriving the correlations is 
exactly the same as before except that now there are four basic correla- 


tions, corresponding to 73; , , components. Without 
dominance, 


“Tee 
é iq aul” 
ag 
‘gy 


STOCHASTIC MATRICES 357 


For example, the transition matrix of full sisters being S;; = 3/33 + 
1T,; , the correlation between full sisters with respect to sex-linked 
traits is 


Tsisters — } + 2(3) i. 


Another example is that the correlation for maternal uncle-niece, with 
transition matrix + 30.; (Table 1), will be 3r;,, = The 
correlations arrived at this way are also in complete agreement with 
those obtained by classical method or by path coefficients. 

14. For sex-linked traits with dominance, the phenotypic correla- 
tions for the four T-component populations are 


Thus, the phenotypic correlation between 


1 (2) 
Full sisters: = 3 + 5 


1 q 
Maternal uncle-niece: r* = 41 


Note that for sex-linked genes, whether with dominance or not, 


Dominance, however, does not effect the correlation between two male 
relatives. 


MULTIPLE ALLELES 


15. The foregoing methods and results may be applied directly to 
multiple alleles, autosomal or sex-linked, without any change except 
that the matrices assume dimensions of 6 X6, 6X 3 or 3 X 3 for the 
case of three alleles with six genotypes. We shall give just one example 
to show the general validity of the method. In the following, the 
order of rows and columns are A,A, , A,;A,,A,A;, A2A2, A2A3, 
for autosomal genotypes or females if the genes are sex-linked. In the 
latter case, the order of males are A, , A, , A; . The gene frequencies 
are p, q, 7, respectively, where p + g + r = 1. The autosomal parent- 
offspring transition matrix is 7,4, = TJ. For sex-linked genes, the 
mother-son matrix is 7; and the father-daughter matrix is T;,. The 
paternal grandmother-granddaughter transition matrix is thus, 


me 
rt = rt = rt =0 
Tas + ’ Tas 1 + Tas 
q q 
| 
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1 0 0 
1 
in gr 0 ® 
7 | 0 P 0 q 0 
03 3 
O 1 
(» q r 0 0 
2p 
29 2p + 2 24 
0 p 0 q r 0 
0 3p + Gr 
LO 0 Pp 0 q r 
We see that exactly the same relation as before holds. It may be also 
verified that Ti; = 47's + 30s. where every row of 0, consists of 


p’, 2pq, 2pr, q°, 2qr, r°. In brief, all of the foregoing relations between 
the matrices hold for multiple alleles. 


DISCUSSION 


16. For autosomal genes, full sibs and double first cousins, being 
able to share both genes identical by descent, are called ‘‘bilineal”’ 
relatives (Cotterman, 1941), while all other types of relatives considered 
here are ‘“‘unilineal” relatives, who can at most share one identical gene 
in common. For sex-linked genes, however, due to the haploid nature 
of the males, the J component appears not only in full sibs of the same 
sex (S;; and S,.) but also in two male relatives with a female as an 
intermediate connecting individual, such as maternal grandfather- 
grandson and maternal uncle-nephew relations (Table 1). 

17. In formulating the relationship between two individuals, how- 
ever, care must be taken to see that the formulation is valid biologically. 
To illustrate what we mean by this, let us consider the autosomal 
relationship given in Fig. 2. We see that the two individuals XY and Y 
are half-sibs and therefore, when X is given, the conditional probabilities 
for Y is 7°. Now, Y and Z are also half-sibs and hence, for specified 
Y, the probabilities of the genotypes of Z are elements of 7°. But, 
when X is given, the probabilities for Z’s genotype is not T°T” = 7", 
as might have been first supposed. Obviously, X and Z are independent 
individuals, for X is the child of A and B while Z is the child of C and 
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D, where the four parents involved are independent individuals. As 
a matter of fact, X is not even related with C, still less with a child of 
C by a second marriage. Therefore, in using the transition matrices, 
it should be borne in mind that one’s relative’s relatives are not neces- 
sarily still relatives. It seems that the rules for using 7 should be the 
same as that for using path coefficients. The latter method is beyond 
the scope of the present discussion; it has been given by Wright 
(e.g. 1934) and explained by Li (1954) in detail. 

18. It should also be pointed out that the subdivision of the correla- 
tion coefficient between relatives as given by r = c; + crrrz is different 
from the subdivision facilitated by the method of path coefficients. 
The latter isolates the components of correlation due to each independent 
common ancestor of the two relatives concerned. Thus, the correlation 
between double first cousins, analyzed by path coefficients, will resolve 
into 4 equal components of magnitude 1/16 each, corresponding tc the 
equal contributions of each of the four grandparents. In this connec- 
tion, particular attention should be directed to the components of 
correlation between full sibs. Either method yields r = } + 3} = }. 
According to our subdivision, the first } is due to identical genotypes 
while the second } is due to half-identical genotypes. By method of 
path coefficients, the first } is due to the influence of one parent (mother, 
say) while the second } is due to the influence of the other parent 
(father). In this special case it is incidental that the two systems of 
subdivision should yield the same apparent results. 

19. One more application of the method of labeling and tracing the 
identity of independent genes may be mentioned. It concerns with the 
genotypic proportions of the offspring of consanguineous matings in a 
random mating population. The half-sibs of Fig. 2 may be employed 
as an example. Let the parents A, B, C be (1, 2), (3, 4), (5, 6), re- 
spectively. The child X may thus be (1, 3), (1, 4), (2, 3), (2, 4), while 


A B c D 


x Y 


FIG; 2. THE INDEPENDENCE BETWEEN THE INDIVIDUAL X AND HIS HALF- 
SIB’S (¥’s) HALF-SIB (Z) 
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his half-sib Y may be (3, 5), (3, 6), (4, 5), (4, 6). If the half-sibs X 
and Y mate, the sixteen possibilities of their offspring may be represented 
by the random union of 1, 2, 3, 4 from X with 3, 4, 5,6 from Y. There 
will be 7 (3, 3) and 7g (4, 4) among the offspring. This 1/8 offspring is 
certainly homozygous because their two genes are derived from the 
same gene of their common grandparent. The probability of its being 
AA is p and being aa is g. The remaining 7/8 offspring is formed of 
independent genes, such as (1, 3), (1, 4), ete. Hence, the offspring of 
half-sib matings consists of 


AA Aa aa 
in agreement with results obtained through other methods. The 1/8 
in the above expression is usually designated by F, known as the 
“<nbreeding coefficient” of the resultant offspring. It is the probability 


that the two genes received by an offspring be identical (derived from 
the same gene of a common ancestor of his parents). 
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A STOCHASTIC MODEL FOR THE SELECTION OF 
MACRONUCLEAR UNITS IN PARAMECIUM GROWTH* 


A. W. KimBa.u anp A. S. HovseEHOLDER 
Oak Ridge National Laboratory 


1. Introduction 


Experimental evidence obtained by various research workers [see 
Sonneborn, (1947) and (1949)] has provided interesting information 
about the existence and behavior of macronuclear units in Paramecium. 
All paramecia of the same stock probably have about the same number 
of macronuclear units. Before division takes place each macronuclear 
unit doubles, and each daughter cell presumably selects one-half of 
the macronuclear units present at the time of division. In this way 
the number of units per animal is kept constant in a manner analogous 
to the behavior of chromosomes. Unlike chromosomes, however, the 
macronuclear units do not separate into homologous sets and as yet 
very little is known about the nature of the macronuclear unit selection 
process. The purpose of this paper is to test one hypothesis about the 
selection process by incorporating it into a stochastic model and by 
comparing predicted results with experimental data. 

In section 2, a state is defined and the transition probabilities are 
derived. Section 3 deals with the theory which underlies the algebraic 
treatment of finite Markov chains. Numerical results for relatively 
small numbers of macronuclear units are given in section 4, and these 
results are compared with experimental data in section 5. 


2. The transition probabilities 


The state of an animal is defined to be the type distribution of its 
macronuclear units (hereafter referred to simply as units). If an 
animal has n units each unit may be different from every other unit. 
or all units may be alike, and any combination between these two 
extremes is also possible. The rth state, defined as a particular distri- 
bution of unit types, may be represented conveniently by the vector, 


(1) a, = , * Ayn) rs om, 


where mis the number of different states and where each a,, 
(i = 1, --- , m) is the number of units of one type. Clearly each ¢ 


*Presented at the American Statistical Association Annual Meeting in Washington, 1953 in a 
session held jointly with The Biometric Society (EN AR) and the Institute of Mathematical Statistics. 
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may assume any integral value from 0 to n subject to the condition that 


n 
a,, =n. 
t=1 


In (1), the order of types is disregarded. For example if n = 3, the 
state (2, 1, 0) implies that the macronucleus contains two units of one 
type, one of a second type and none of the third type. No distinction 
is made between two macronuclei, one having two A-units, one B-unit 
and no C-units and the other having two B-units, one A-unit and no 
C-units. Two animals having such macronuclei are said to be in the 
same state, (2, 1, 0). 

Readers familiar with number theory will recognize that the number 
of possible states, m, is identically equal to the number of partitions of n. 
Thus the number of possible states is the coefficient of ¢* in the power 
series expansion of 


T]a- 


No workable formula for m as a function of n has ever been derived 
although much work has been done on asymptotic formula and re- 
currence relations. 

For certain experimental problems information is needed as to the 
change in the distribution of states in a Paramecium population as the 
population grows. If a population is generated from a single animal 
whose state is known and if the probability of transition from one state 
to any other state as a result of division can be computed, the distribu- 
tion of states in the population after v divisions can be obtained by 
treating the growth process as a finite Markov chain. Likewise if a 
population is generated from a group of animals and if the initial distri- 
bution of states is known, the distribution after v divisions may be 
derived directly from the results obtained for a population generated 
from a single animal. The development throughout the paper assumes 
that all animals in the population undergo the same number of divisions 
per unit length of time. While this is most certainly not true for any 
specific time period, experimental evidence indicates that the assump- 
tion becomes more closely approximated as the period of observation 
increases. In principle division rates could be incorporated into the 
problem by approximating the distribution from experimental data and 
by using this information in establishing the transition probabilities. 
Whether this approach would be feasible is questionable. 

If an animal normally has n units in its macronucleus, then at the 
time of division 2n units are available to the daughter cells each of 
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which ultimately receives n of them. One possible mechanism for the 
manner in which a daughter selects n out of the 2n units is the hypo- 
thesis of completely random selection. According to this hypothesis 
each of the (°*") ways of selecting n units out of 2n units is assumed to 
be equally probable. For the purpose of computing transition probabili- 
ties only one daughter need be considered, since the second daughter’s 
lot is uniquely determined after the first daughter’s selection is made. 
In general, the probability of transition from state a, to state a, is 
given by 


(2a,,)! (2a,s)! +++ (2a,,)! 
n 


where a, 8, --- , o are the numbers of like a,; and 


1 
The symbol | z,,; | refers to the permanent, a mathematical form 
identical to the determinant except that all terms are taken with the 
positive sign. In (3), if a,; > 2a,; , x;; is taken to be zero, consistent 
with the usual definition of the factorial of a negative integer. 

The validity of (2) can be proved rigorously but it can also be 
deduced almost directly from a simple example. Consider the case 
n = 3. Here m = 8 and the possible states are (1, 1, 1), (2, 1, 0) and 
(3, 0, 0). In general we are interested in the probability of transition 
from state a, = (a@,; , @,2 , @,3) to state a, = (a,, , G.2 , G3). If the 
elements of a, are all different, the number of ways state a, may be 
obtained from state a, is 


Qa As2 Qs3 Asi aay 
(4) 
as; As3 as; Aso 53 


which may be written in the form, 


where the z,; are defined as in (3). The expression in brackets ii: (4’) 
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+ 
is immediately identified as the permanent | x,,; | , since it contains 


exactly the six terms of the corresponding third order determinant 
but with all signs positive. The expression in (4’) was derived on the 
assumption that the a,; were all different. If in fact some of them are 
alike, we must make the usual correction which consists of dividing 
by the factorials of the numbers of like a,; . Finally after division by 
(), We arrive at the desired probability given in (2). The extension 
to cases in which n > 3 is almost intuitively obvious. 

The statement that the transition probabilities can be computed by 
considering only one daughter may not at first glance appear valid. For 
example if n = 4, there are five possible states, (1, 1, 1, 1), (2, 1, 1, 0), 
(2, 2, 0, 0), (3, 1, 0, 0) and (4, 0, 0, 0). An animal in state (3,1, 0, 0) 
can upon division give rise to daughters which are in states (2, 2, 0, 0), 
(3, 1, 0, 0) or (4, 0, 0, 0). If from the six units of type A and 2 units 
of type B which are available one daughter selects 2 of type A and 2 of 
type B, the remaining daughter will be left with 4 of type A and none 
of type B. That is, the two daughters will be in different states, 
(2, 2,0, 0) and (4, 0, 0, 0). This might appear to make our calculations 
invalid. Fortunately, however, each result of this kind is exactly 
balanced by another selection which occurs with exactly the same fre- 
quency. In this case if the first daughter selects 4 of type A, the second 
daughter is left with two of type A and two of type B, so that the two 
daughters are again in the states (2, 2, 0, 0) and (4, 0, 0, 0). The net 
result is that the relative frequency of occurrence of each state is exactly 
what it would be if daughter pairs were always alike. For many tran- 
sitions, of course, the two daughters from a single animal actually are 
in the same state. 


3. The general theory 


Since reproduction does not increase the number of distinct types uf 
units, a daughter cell can have at most as many types as the parent, and 
there is a finite probability that the number will be smaller. The total 
number of possible states is m, the number of partitions of n, as indicated 
already. These states can be divided into n groups, according to the 
number of distinct types of units present. Let the groups be designated 
9: » *** » Jn , Where g; is the group of those possible states for which 
exactly 7 distinct types are present. Then no descendant of a cell of 
group g; can be in a group g; for 7 > 7. 

The vector a, = (a,; , @,2, *** , Gn), Where the a,; are non-negative 
integers satisfying 


=n, 
i 
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is taken to represent the state in which there are a,, elements of one 
type, a,. of another, etc. It is no restriction to suppose that 


We can make a lexicographic ordering of the possible states and let the 
subscripts r correspond to this ordering, as follows. If r < s, then for 
some 17, 


0 = Gin = = °°? = G41 > a,; — 


Thus in state a, , all elements are alike; in state a, , all but one are 
alike; in state a; , there are two elements of one type and n — 2 of 
another; --- ; in state a, , the elements are all different. We are in- 
terested in determining, as a function of v, the probability that a daughter 
cell of the vth generation is in group g, , or, more generally, in g; or 
some lower group. 

Let p represent the matrix of probabilities p,, . Let po;, represent 
the probability that some particular cell is in state a, , and let po) represent 
the row vector whose elements are these probabilities. Then 


Pi = PoP 
is a row vector whose sth element is the probability that a particular 
daughter cell is in state a, .. More generally, 

= Pop” 


is a row vector whose sth element is the probability that a particular 
daughter cell of the vth generation is in state a, . We recall that all 
elements of p and all elements of each p, are non-negative, and that 


Pris = = 1, 


The classification into groups of states provides a natural partitioning 
of the matrix p and, in fact, it has the form, 


p= Qu GQ O 
Q32 Qs P 


where the submatrices above the diagonal are all null because of the 
impossibility of a transition to a higher group. The principal submatrix 
q: is a scalar and, in fact, gq, = 1 since group g, consists of the single 
state in which all elements are alike. The order of q, is either n/2 or 
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else (n — 1)/2, whichever is an integer. Both q, and q,-_, are scalars, 
Qn-2 1s of second order, and q,_; is of third order*. 

If, for any square matrix A, we let A(A) represent the characteristic 
function, then in this notation, 


P(A) = gu(d). 


Since ¢,(A) = 1 — \, it follows that \ = 1 is a characteristic root of p. 
In a matrix of positive elements, the greatest characteristic root cannot 
exceed the greatest sum of the elements of any one row and cannot be 
exceeded by the least of these. But every row of p beyond the first 
has at least one non-null element that is not an element of any q,; . 
Hence, q;(1) ¥ 0 for every i > 1, and 1 is therefore a simple root of 
p(A) = 0 and exceeds in magnitude every other root. Thus if \, = 1, 
are the distinct characteristic roots of p, 


Let A represent the Jordan normal form of p, so that 
p= M~'AM 


for some non-singular M. In A each diagonal element is a \, , the first 


subdiagonal consists of ones and zeros, and all other elements are null. 
Then 


p’ = M"'A’M 
and 
= ppM'A’M. 


Each element of the vector p, can be expressed as a linear combination 
of the A; (¢ = 1, --- , u) whose coefficients are either constants or are 
polynomials in v. If A is a diagonal matrix the coefficients sre constant. 
If X, is a root of multiplicity m, , then the coefficient of A; is a poly- 
nomial in v of degree at most m,. Hence 


= ky + + + , 
where the k,(v) are vectors whose elements are either constants or 
polynomials in v. 

For increasing v each row of p’ approaches a characteristic row 
vector associated with the largest characteristic root of p, and each 
column approaches a characteristic column vector associated with this 
same characteristic root. However the largest characteristic root is 


*Forn < 6,¢,-, and q,_, may be of lower order. 
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, = 1, and the associated characteristic row vector is the unit vector 
= (1, 0, , 0). 


Since p, — k, as v © and since the sum of the elements is always 
unity, it follows that k, = p, , so that 


Pv = Pi + t+ 


We are interested in the rate at which p, — p, . If there is a root 
\» Which exceeds in modulus all the other roots except \, then asymp- 
totically, 
= pi + . 


For the special cases n < 7 it turns out that A, = (2n — 2)/(2n — 1), 
and this is a simple root. Hence for these cases, k,(v) is a constant 
independent of v. This root satisfies g.(A) = 0. For n < 6 it is also 
true that the next largest root is 43 = (2n — 4)/(2n — 1), which is 
again a simple root and satisfies g,(4) = 0. The matrix g, may be 
expected to have the largest root different from unity since its row sums 
will be greatest. It remains to be determined whether these formulas 
are valid generally. 

Although it has not been determined that 4, = (2n — 2)/(2n — 1) 
in general, nevertheless we can set a bound for this root. The matrix 
Yo: is, in fact, a column vector whose jth element can be written 


(?") 2n(2n 1) eee (2n 23 1)’ 
n 


The smallest of these elements is the last for which 7 = [3], and is equal 
to 1/(*) if n is even and (n + 1)/(*) if n is odd. If we call this y, , 
then 1 — y, is the greatest row sum for the matrix g, and 1 — (n — 1)/ 
(4n — 2) = (8n — 1)/(4n — 2) is the least. Hence d, is bounded 
between these two limits.* 

For more complete information one would like, of course, to have 
general formulas for all the \’s and the corresponding k’s. Intuitively 
the nature of the model would seem to rule out the possibility for 
oscillatory behavior, and we might therefore expect that all \, would 
be positive and real. This is borne out in all cases examined (n < 6), 
the roots actually being positive rationals. A few of the roots have been 
determined in general: q, and q,_, are both scalars having the values 


*Dr. Motoo Kimura, one of the referees, has shown that A: = (2x — 2)/(2n — 1) holds asymptoti- 
cally. 
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n n 
Qn-2 is of order 2 and has the roots 
+ ?) 
( 2) 
(?") 
n n 
and q,-3 is a matrix of order 3 and has the roots 
n—1 n + ') + 3) 
( 1 
n n n 
The roots of g; , as 7 decreases by unit steps from 7 = n, form a pattern 
which suggests that the roots of g,-, would be 


and similarly for the remaining g; . Although we can offer no proof of 
this, the formulas do provide the correct roots for n < 6. 


4. Numerical results for some special cases 


To illustrate how the computations are carried out, we shall consider 
the case n = 5 in detail. In this case m = 7, and the possible states are 


a, = (5, 0, 0, 0, 0) 
a, = (4, 1, 0, 0, 0) 
a; = (3, 2, 0, 0, 0) 
a, = (3,1, 1,0, 0) 
a; = (2, 2, 1, 0, 0) 
a, = (2, 1, 1, 1, 0) 
a, = (1,1, 1,1, 1). 
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The first step is the calculation of p, the matrix of transition probabili- 
ties. From (2), for example, we find 


11111 
11111 
10 212121212! 
= 1111 1) = 82, 
and 
1 1/6 1/6 1/24 1/24 
01 1 #12 1/72 
10 412121210! 
= O1 1 12 1/2 | = 48. 
01 1 12 12 
00 0 1 1 


In evaluating the permanents, two properties are useful. Constants 
may be factored from rows or columns as in determinants, and the 
permanent of an n X n array consisting of all ones is n!. Repeated 
application of (2) to obtain all the transition probabilities yields 


(252 0 0 0 0 0 0 


56 140 56 O 0 0 0 
66 180 0 0 0 0 
1 
P = 559 60 40 80 66 0 O}. 
0 
0 


6 12 48 90 9% 
0 0 0O 60 160 32) 


The next step is the determination of the characteristic roots and 
vectors of p. Three roots may be obtained directly and the other four 
are solutions of the two equations, 


6 
6 
0 12 56 64 120 
0 
0 


(80/252 — d) 66/252 
64/252 (120/252 — 2) 

| (140/252 — 2) 56/252 
66/252 (180/252 — 2) 
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Hence, the roots* are 
\Ar Ao, Ax) = 1/252(252, 96, 224, 32, 168, 96, 32). 


The characteristic vectors, which may be obtained in numerous ways, 
are 


[‘m,| 1 0 0 0 0] 
Mm: 1-3 2 0 0 
ms -2 ll 14 0 oOo 
M = |m| = of o-b 4i, 
ms 19 -27 -30 16 22 0 0 
Me 0 oO 1 <1 1 0 
im) 5 -10 0 30 -40 16) 


where m, is the characteristic vector corresponding to \,; . 
At this point, given any initial probability distribution of states 
Po , We could compute the distribution after v generations, 


p, = pM"a’M, 


after first computing the inverse of M. (In this case A is a diagonal 
matrix the elements of which are the A; .) For the purpose at hand 
(see section 5), we are interested only in the special case in which all 
cells initially are in the state a; , i.e., the state in which all units are 
different. This corresponds to the initial distribution pp = (0, 0, 0, 0, 
0,0, 1). Thus it is clear that we need only the last row of the inverse 
of M which is 


pM = (1, —35/64, 5/64, 9/68, 15/204, 5/2, 1/34). 


Furthermore (see section 5), we are interested only in the first element 
of p, which is the probability corresponding to the state a, . In other 
words this element gives the probability that, after » generations, a 
daughter cell which originally came from a cell in which all units were 
different would have all its units alike. The element is 


1 — (35/64) (96/252)” — (125/64) (224/252)’ + (9/68) (32/252)’ 
+ (285/204)(168/252)’ — (1/34)(32/252)’. 


Expressions giving this probability have been derived for n = 2, 


*Ordered to correspond with the rows of p. 
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3, 4, 5 and are shown in Table 1. The corresponding survival curves 
have been computed and are shown in Figure 1. The dotted lines in 
Figure 1 represent experimental data which are discussed in the next 
section. 


TABLE 1. 


The probability that, after v generations, a daughter cell will have all units alike, given that the original 


parent cell had 211 units different. 


Units per Probability 


Cell (n) 


— (2/8)’ 

(3/2)(4/5)” + (1/2)(2/5)’ 

(39/22) (30/35)” + (20/35)” — (5/22)(8/35)" 
(125/64) (56/63)" + (95/68) (42/63)” 

— (35/64)(24/63)" + (7/68)(8/63)" 


| 


5. A comparison with experimental data 


The data which provide a test of the proposed model for macronuclear 
unit selection actually came from an experiment with Tillina carried 
out* by Dr. Josephine Bridgman. It seems certain that the macro- 
nucleus of Tillina is of basically the same structure as that of its near 
relative Paramecium. ‘Thus we are justified in using the data from this 
organism for our comparison. 

In the experiments considered here, tillinas were irradiated at a dose 
which resulted in essentially 100 per cent mortality after about twelve 
divisions. Such severe early mortality suggests that each of the n 
units was rendered lethal as a result of irradiation. If this is true and 
if we assume that all animals started out in the state (1, 1, -*- , 1), the 
survival curves which would be expected for various numbers of units 
are precisely those shown in Figure 1. That is, these curves represent 
the expected proportion of a population which would not be in state 
(n, 0, --- ,0) after divisions. If each unit is lethal the state (n, 0 ---,0) 
would be synonomous with death. 

The results of two experiments plotted as per cent survival against 
number of divisions are shown as the dotted lines in Figure 1. If the 
model and assumptions are correct, the number of macronuclear units 
would have to be in the neighborhood of 3 or 4. Since the actual number 
of units is likely to be 20 or greater, the model or the additional as- 


*Work performed in the Biology Division, Oak Ridge National Laboratory. 
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sumptions, or both, must be discarded. Clearly any change in the 
assumption concerning the initial distribution of states would result in 
a higher estimate of n. That is, if any of the original animals already 
had some of their units alike, the predicted survival curves would be 
uniformly lower since in each case the state (n, 0, --- , 0) would be 
reached more rapidly by the animals starting out with some units alike. 
This would be one way to obtain a better fit to the model. Actually, 
however, among less than one hundred units, it is not very likely that 
any two would mutate lethally in exactly the same way. For this 
reason the assumption that all animals are initially in the state 
(1, 1, --- , 1) is probably a close approximation to reality. 

One reasonable assuraption which would bring the theory and 
experiment closer together is that some number of like units less than n 
would produce death. : Previous experience with other organisms has 
suggested that with recessive lethal mutations and a high degree of 
polyploidy, death can result from states of less than complete homo- 
zygosis. Thus, if we assume that death can result if all but k of the 
nm units are alike, k may be varied until the theoretical curve for, say, 
n = 20 approaches the experimental curve. If the value of k required 
to bring the two curves together is too large, say greater than 50 per 
cent of n, this approach would also have to be discarded since it is not 
likely that death would result if only half of the units were alike. These 
computations will not be practicable until more complete results are 
obtained for the general case. 

Finally, one could consider various changes in the basic model. 
Allowance might be made for spontaneous mutations, or some hypo- 
thesis other than complete random selection might be proposed. Un- 
fortunately, in this experiment some evidence is available which suggests 
that the macronuclear unit selection hypothesis may not be the proper 
way to describe the lethal effects of radiation. In some respects the 
death phenomenon seems to be an all-or-none affair in the sense that 
death for an animal in any particular line of descent can be predicted 
rather well after the first division has occurred. Furthermore, division 
rates are not uniform and seem to depend on whether the animal is 
destined to expire. If these observations are correct, the proposed 
model may never fit the data very well, and even if it does, its validity 
would be questionable. Nevertheless, the stochastic process itself is 
rather general in nature and perhaps has counterparts in other biological 
problems or even in other fields of endeavor. 
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6. Summary and conclusions 


A stochastic model is proposed for the selection of macronuclear 
units in Paramecium growth. The model presupposes random selection 
and that death can result only in the completely homozygotic state. 
Although full results for the general case have not been obtained, 
enough has been done to indicate that the model must be changed if 
reasonable agreement with experimental data is to be achieved. The 
most plausible alteration which would provide better agreement is an 
allowance for death in states somewhat less than complete homozygosis. 
Before computations to test this hypothesis can be started, more 
results are needed for the general case of m macronuclear units. At 
present complete results are available only for values of n up to 6. 
Hand computation for larger values of n rapidly becomes impracticable. 

The authors are indebted to Dr. R. F. Kimball, Biology Division, 
Oak Ridge National Laboratory, for suggesting this problem and 
assisting in its formulation and for a careful review of the manuscript. 
We also wish to thank Dr. Bridgman for permission to use her hitherto 
unpublished data. 
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INCOMPLETE BLOCK RANK ANALYSIS: ON THE 
APPROPRIATENESS OF THE MODEL FOR A 
METHOD OF PAIRED COMPARISONS.’ 


ALLAN BRADLEY 
Virginia Agricultural Experiment Station of the 


Virginia Polytechnic Institute 


1. Introduction: 


The application of any statistical technique to numerical data is 
strictly appropriate only when the data conform to the assumptions 
and requirements inherent in the development of that statistical 
method. In the analysis of variance, we require normality, additivity, 
independence, and homogeneity of variances of the data. When these 
conditions are not at least approximately satisfied, the analysis is 
suspect. In postulating the mathematical model for the rank analysis 
of incomplete block designs and more specifically for the method of 
paired comparisons [1] we chose the model because it seemed intuitively 
plausible, because it was related to the binomial model, and because 
it was mathematically workable. 

The major objective of this paper is to develop a procedure for 
testing the appropriateness of the model for the method of paired 
comparisons. The test is then applied to a variety of experiments 
involving taste, preference, and appearance judgments for a technique 
is only very useful when it is applicable to a considerable variety of data. 

In the reference cited above, three tests of significance were proposed. 
Two of them are tests of treatment effects, and the other is a test of 
agreement among judges or over groups of the data. The test of agree- 
ment could be regarded as a test of the need of the model that permits 
judges to differ on their judgments as opposed to the model in which 
it is postulated that all judges are coordinated in their judgments on 
the attributes under consideration. We could call the test of agreement 
a between-group test of goodness of fit. 

These comments will become clearer in view of the following 
summary of the mathematical model for paired comparisons. 


'Prepared under a Research and Marketing \et Contract, Project RM: e-629-1, with the Bureau 
of Agricultural Reonomies. 
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2. The model and tests of hypotheses: 


Let us consider the method of paired comparisons. A basic repe- 
tition of a paired comparisons experiment may be defined to be a set of 
incomplete blocks of size two in which every pair of ¢ items are com- 
pared once and only once. The experiment is taken to consist of n 
such basic repetitions. Various procedures for the analysis of paired 
comparisons are available. The method of analysis used depends on 
the form in which data are recorded: measurements or scores may be 
available for each sample or for the difference between items in each 
pair, or only the preferred item in each pair may be indicated, thus 
constituting a ranking procedure. Other innovations, such as the 
grouping of repetitions of the experiment by judges or by time or by 
some other relevant characteristic, may be developed for further experi- 
mental control. We shall review the formulation of the mathematical 
model and the test procedures for a method of paired comparisons, the 
rank analysis of incomplete block designs, wherein items are only 
ranked and some control grouping of repetitions of the experiment may 
be present. 

In summarizing the model for the method of paired comparisons 
under consideration, we shall first assume that the experiment consists 
of n repetitions for ¢ treatments and that the repetitions are a homo- 
geneous set. Then it is postulated that associated with each of the ¢ 
treatments, say T, , --- , JT, , there exist parameters, z; for T; , such 
that +; = 0 and >°‘_, x; = 1. The latter restriction is one that is 
imposed to insure determinancy. The behavior of the parameters is 
further defined with the probability statement that, if X; generally 
denotes an observation on a sample of 7; , 


(1) P(X; > Xi) = + 17) 


in the comparison of 7; with T; . 

Experimental observation is limited to ranking in the pairs. Thus, 
r,;x 18 defined to be the rank of 7; in the block in which 7’; is compared 
with 7; in the kth repetition of the experiment, 1 S k S n. A non-null 
decision is required in each pair and of necessity r;;, is either one or 
two with r;;, + 7; = 3. The rank of one will normally be assigned 
to that treatment of a pair that is judged ‘superior’ on the basis of the 
test attribute of the treatments. 

The element of the likelihood function for the block containing 7’; 
and 7’; may be written as’ 
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The general likelihood function for the ¢ treatments and n homogeneous 
repetitions is abbreviated through use of the functional form L(x,) with 

where the index of summation or of multiplication covers the full range 
when not otherwise specified, where the prime indicates that the index 
of summation j takes all values except 7, and where the second product 
contains one factor representing each treatment pair. The maximum 
likelihood estimator of 7, is denoted by p; . Procedures for obtaining 
these estimators are outlined in two papers [1, 4]. We shall need to 
refer to (2) in developing the theory for testing the appropriateness of 
the model. 

In some cases the model that assumes homogeneous repetitions 
is not realistic. This is particularly true in certain types of taste- 
testing experiments where a number of judges perform repetitions of the 
experiment and where the main objective is the detection of treatment 
differences. Then it may be reasonable to postulate the existence of 
parameters , , Tu, Tix = 0, = 1 for the uth of g groups 
of repetitions in the experiment. This simply means that the experi- 
ment may be designed to contain g groups of repetitions with n, repe- 
titions in the uth group, >>, », = n. Repetitions within a group are 
taken to be homogeneous, but group differences in treatment para- 
meters are permitted. The mathematical model differs from the simpler 
case described only in the grouping of the repetitions. The likelihood 
function in this case is a simple product of g functions of the form (2). 
The maximum likelihood estimator of 7;,, may be designated by pj. . 

The test procedures that have been developed under these models 
for paired comparisons are as follows: 


Test (i): H,: «; = 1/t for all z, 
H,: # 1/t for some 7. 
The likelihood ratio test depends on the statistic? 
(3) log + pi) — {2n(t Tis log Ds 
When we let \, be the likelihood ratio statistic for test (2), 
(4) —2 Ind, = nt(t — 1) In2 — 2B, In 10 


has for large values of n the x°-distribution with (¢ — 1) degrees of 


?We shall use log and In to represent logarithms to base 10 and base e respectively. 
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freedom*. This first test is a test of treatment equality under the 
assumption of homogeneous repetitions. 


Test (ii): for all i and u, 
for some 7 and u. 


We define B! to correspond to B, of (3) for the n, homogeneous repeti- 
tions of the wth group and 


(5) Bi = 
The likelihood ratio statistic \; for Test (ii) is defined by 
(6) —2 Ind; = nt(t — 1) In2 — 2Bi In 10 


which has for large values of n, , --- , n, the x’-distribution with 
g(t — 1) degrees of freedom‘. 


Test (iii): for all 7 and u, 
for some 7 and wu. 


Test (ili) may be regarded as a test of agreement of treatment para- 
meters over the g groups, and this is analogous to a test of treatment X 
group interaction. In this case, we may define the likelihood ratio 
statistic \, by 

(7) —2 Ind, = —2(In — IndA,) 


2(B, — Bi) In 10 


where A, is obtained by assuming that all g groups of repetitions may 
be pooled to form an experiment with n homogeneous repetitions. The 
large-sample distribution of —2 In \, is that of x’ with (g — 1) (¢ — 1) 
degrees of freedom, and the distribution for small samples is dependent 
On , , aS nuisance parameters. 


3. Tests of goodness of fit: 


We shall limit our main consideration to a test of goodness of fit 
for the simple paired comparisons model for ¢ treatments and n homo- 
geneous repetitions of the possible pairwise treatment comparisons. 


3The exact distribution of 4, has been tabled for small values of and xn. Tables fort = 3, = 1, 
*** , 10 and fort = 4,2 = 1,°** , 6 have been published [1]. Unpublished tables for ¢ = 4, n = 7, 
8, and fort = 5,n = 1,°** , 5 have been prepared. 

‘Limited tables for the distribution of B,° are published [1]. Only cases where mi, *** , Ng are all 


equal are considered, 
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It will then be easy to extend this theory to the more general model 
with grouped repetitions. 

The most general model available for a paired comparisons experi- 
ment within a group of data is one in which each comparison is treated 
independently as a small binomial-type experiment. One parameter 
is then estimated for each comparison. For example, when treatments 
T; and 7’; are compared, a parameter 7,; , the probability that 7’; is 
ranked above 7’; , would be estimated. The complementary probability 
;; = 1 —7,; . Now, if ¢ treatments are involved in an experiment of 
this type, (:) parameters must be estimated. Estimates are obtained 
from relative frequencies so that 


(8) Dis = fis/n 

where p;; is the estimator of 7,;; , f;; is the number of times that T; 
ranks above T; , and n is the number of times such a comparison is 
made. This general model fits the experimental data exactly in the 
sense that observed and expected frequencies of a particular ranking 
are identical.° When (;) parameters are to be estimated, the likelihood 
function may be written as 


where 
fi and + = 1. 
To test the appropriateness or the goodness of fit of the model for 


paired comparisons for homogeneous repetitions, we are interested in 
the following test: 


Test of Fit I: Hy: = + 7;), sit 
H,: + for some and ). 


An approximate test is possible using the likelihood ratio procedure. 
The likelihood ratio statistic for Test of Fit I depends on L(j;; | Ho) 
and L(p;; | H;) where L is defined by (10) and these functions represent 
evaluations of L in terms of estimators p;; obtained under the assump- 
tions of H, and H, respectively. Further, it will be noted that 
L(pi; | Ho) = L(p,) where L is defined in (2). These quantities are 


‘The more restrictive model of Section 2 is obtained when we set #;; = #;/(*i + #;). Further, 
the relative frequencies are related to the sums of ranks previously detined for 


When it can be done without ambiguity, we shall use © r; as an abbreviation for the left-hand member 
of (9) and refer to it as the ‘total sum of ranks for 7;’. 
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evaluated using the appropriate maximum likelihood estimators and 
we have 


(11) In L(p;; | H.) = —B, In 10, 
and 
(12) In | H;) = fii In fii n( £) Inn 


in view of (3) and (8). B, is the statistic tabled for Test (i) and is 
available without special calculation for a fairly wide range of values 
of ¢and of n. For the test of goodness of fit, the statistic, 


p> n( Inn + B, In 10, 
iti 

has the x’-distribution with (;) — ¢ + 1 degrees of freedom for large 

values of n. 

The computation for the test is most easily carried out after forming 
the two-way table commonly used in tabulating the results for methods 
of paired comparisons. For three items or treatments, this table appears 
as follows. 


TABLE I. TABLE FOR THE TEST OF GOODNESS OF FIT 


T: Ts 
Siz his Zn = 4n — fix — fis 
- Sas = 4n — fu — fos 
Ts Su Su = 4n — fu — fn 


When the repetitions of the paired comparisons experiment are 
grouped into g groups, the uth of which contains n, repetitions, uv = 1, 

-,9, Luu % = Nn, separate and possibly distinct sets of parameters, 
7.;, and 7;, are assumed to exist within each group. Statistics f;;, 
and >> r;, are computed for each of the g groups. Within-group tests 
of goodness of fit are made as described for Test I and the test statistics 
may be designated by —2 In We note that all of these parameters 
and statistics correspond to those of the simpler case and that we have 
only added the additional subscript u to designate the groups of the 
experiment. 

The test specification for an over-all test of goodness of fit is 
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Test of Fit IT: 
Hay: A + for some i and j and u. 


In view of Test of Fit I and the fact that rankings by repetitions of the 
experiment are taken to be independent in probability, the likelihood 
ratio statistic for Test of Fit IT is 


(14) = —2 >> Ind, , 

and this statistic has the x*-distribution with g{(3) — ¢ + 1} degrees 
of freedom for large values of 7, . 

Test of Fit IT should be used when repetitions of the experiment 
are grouped. When grouping is employed, it must have been necessary 
to assume a priori that one set of parameters, 7, , --- , 7, , would not 
be sufficient to characterize the experimental results. It would not thus 
be reasonable to use Test of Fit I and to ignore the grouping in these 
situations. 


4. The tests of goodness of fit associated with expected frequencies: 


In tests of goodness of fit, the procedures usually involve computing 
expected cell frequencies, say f!; , and forming x° by taking sums of 
terms of the form, (f;; — f/;)° f/; .. Again, in this case, we can show 
that our tests, which depend on logarithms of cell frequencies, are 
approximately equivalent to the more usual ones. We shall consider 
Test of Fit I, and the results obtained will apply directly to the second 
test as well. 

Sums of ranks are directly related to observed cell frequencies 
through (9) and Table I. Expected cell frequencies are related to the 
estimators p, , --- , p, through 


(15) fi; = nps/(pi + Pi). 
Then, with the use of (3) and (13), B, = — ) a fi; log (fi, /n), and 
(16) —2Ind; = 2 fi; In 

It will be convenient for the moment to write f;; f{; = 1+ «; , 


where ¢,; may be either positive or negative. Now, 


=2 +) mn (1 + 
We use power series expansions for the logarithms stopping with the 
second terms. The errors committed in so doing will not be large 
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if | «;; | is small. Difficulty is introduced, however, if €;; = — 1 which 
occurs when an observed cell frequency is zero. Upon expanding the 
logarithms and noting that >-,.; f/; «;; = 0, we have 

—2 In => fie; — fini: 


1*) 


The final result is that 


if we neglect the terms involving ¢:; and note the relationship of ¢,;; 
with cell frequencies. 

The procedure for the second test of goodness of fit is clear since 
the test statistic is made up of a sum of terms like the left-hand member 
of (17). We shall retain again the definition of the symbols used in (17) 
but add another subscript u to specify the group of repetitions. Then. 


(18) —2 In Ay = > 


The two methods of computing x’ for tests of goodness of fit or of 
the appropriateness of the models will usually be equivalent for practical 
purposes. Further, there is little difference in computing time required 
when values of B, are available. Otherwise, the second forms, (17) and 
(18), may be slightly easier to handle. Both methods require values of 
Pi, °*** , Pe, Which must be obtained from tables or computed for the 
experiment. The additional work in making these tests of the basic 
models over the analysis of the data for test purposes is very small. 

We have noted that difficulty is encountered in obtaining the forms 
(17) and (18) if observed cell frequencies are small. In addition, the 
x’-tests developed are only approximate and require reasonably large 
numbers of repetitions. If one were about to appraise the applicability 
of the model of the method of paired comparisons to a specific type of 
experimental problem on which a considerable amount of experimenta- 
tion is planned, it would be very useful to conduct a fairly extensive 
experiment to obtain data for tests of goodness of fit. 


5. Applications of the tests of goodness of fit to experimental data: 


In the following subsections, we indicate the results of testing the 
appropriateness of the model postulated for the rank analysis of in- 
complete block designs for three different types of data. These sets of 
data are those that happen to be available in the form required for the 
test. No deliberate selection has been made, and the reader may 
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assess the results as they were obtained. The data were obtained by 
the credited research workers exactly as shown, and we have only 
obscured the quantitative nature of the treatments because the data 
have not yet been published elsewhere. 

It does seem necessary that some treatment differences exist before 
data can show any marked departures from the model used. Thus it 
is to show that these are not simply ‘blank’ experiments that the levels 
of significance for tests of treatment effects have been included under 
each of the following tables. These significance levels should not be 
confused with significance levels associated with values of x’ also given 
under the tables. 

(i) Experiment 1. The Effects of Peanuts in Hogs’ Rations on 
the Flavor of Fresh Pork Roasts. C. M. Kincaid, H. R. Thomas, and 
Lyle L. Davis, Virginia Agricultural Experiment Station. 

We shall present the data from this experiment summarized in 
tables similar to Table I. The treatments consist of different rations 
fed to hogs that were selected for similar characteristics and from a 
limited number of litters. A is a corn ration, B is a corn ration supple- 
mented by a small peanut ration, and C is a corn ration supplemented 
by a larger peanut ration than was the case for B. The experiment 
was conducted so that the tasters recorded a rank of 1 for that sample 
in a pair that they thought came from the animal fed the larger amount 
of peanuts. Members of the panel were obtained after preliminary 
training and selection, and it was believed that they would be able to 
recognize the effect of a peanut supplemented diet on fresh roast pork. 
Actually these data do not reveal any over-all treatment effects but are 
somewhat indicative of detectable differences among roasts in the sets. 

During the conduct of the experiment, even after the preliminary 
training, it appeared that there was still some difficulty in coordinating 
the rankings of the judges. For this reason, the combined analysis 
[Test (ii)] was used for tests of treatment differences and we now test 
the model for each judge separately. It also appeared that differences 
among roasts might be as important as differences between rations, 
and hence we also test the model for each set of roasts over the panel of 
judges. Five judges and five sets of roasts were used in the experiment. 

The required tables for these ten tests of the model are shown below. 
The calculations required are illustrated in detail for the first judge, 
and we note that each x? computed in Experiment 1 has 1 degree of 
freedom. In each case here x’ has been calculated using both of the 
forms (13) and (17) with the latter value being shown in parentheses. 
The 5 per cent significance level of x? with 1 degree of freedom is 3.84 
and the 1 per cent level is 6.63. 
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TABLE II. RESULTS FOR JUDGE I 


A B Dr; 
A 7 5 28 
(6.00) (6.00) 
B 3 ~ 6 31 
(4.00) (5.00) 
5 4 31 
(4.00) (5.00) 


t = 3,n = 10, B, = 8.8560, x? = 1.24, (1.23) 
Sig. level of treat. = .78. 


Using Test of Fit I in the form (13), we have 


—2 Ind; = xi = 2(2.3026)[7 log 7 + 5 log 5 + 3 log3 
+ 6 log 6 + 5 log 5 + 4 log 4 — 30 log 10 + 8.8560] = 1.24. 


This value is shown below Table II followed by the value obtained 
using expected cell frequencies. 

To use the second form (17), we require p, = .4286, pg = .2857, 
and pc = .2857 obtained from tables [1]. The expected cell frequencies 
are shown in parentheses in Table II and they were calculated using 


the relationship, (15). The alternative method of computing x’ yields 


(7-6) , 3 — 4)’ 
+ 
6 6 4 
(6— 5)? , 6-4) , 
= 1.23. 
+- 5 + 4 + 5 
TABLE III. TABLE IV. 
RESULTS FOR JUDGE II RESULTS FOR JUDGE III 
A B C A B Dr; 
A - 3 4 33 A - 5 3 32 
4. B 7 - 5 28 B 5 - 5 30 
6 5 29 7 5 28 
t = 3,n = 10, B, = 8.6191, t = 3,n = 10, B, = 8.7973, 
xz = 0.15 (0.16) x? = 0.57 (0.57) 


Sig. level of treat. 


50. Sig. level of treat. = .63. 
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In all cases considered, the two methods of computing x’ for Test of 


Vit I yielded values in close agreement. 


TABLE V. 
RESULTS FOR JUDGE IV 


ari 

A 3 33 

B 7 —- 5 28 
5 


C 6 | | 29 


t=3,n= 10, B, = 
x? = 0.15 (0.16) 
Sig. level of treat. = .50. 


| 


TABLE VI. 
RESULTS FOR JUDGE V 

= | 
A B Cc 

| 
A _ | 4 4 32 
B 6 _ 8 26 
C 6 2 ~ 32 


Sig. level of treat. 


t = 3,n = 10, B, = 8.3162, 


xi = 1.37 (1.35) 
.25. 


ll 


Test of Fit II may be applied here to the over-all experiment. 


Then, using (14), we have 


—2IndAy = x3 = 1.24 + 0.15 + 0.57 + 0.15 + 1.37 = 3.48 


TABLE VII. 
RESULTS FOR ROASTS Ai, Bi, Ci 


A B C zr; 


TABLE VIII. 
RESULTS FOR ROASTS Bz, Cz 


A — 5 3 32 
B 5 = 2 33 
C 8 25 


t = 3,n = 10, B, = 7.8787 
xi = 0.19 (0.19) 
Sig. level of treat. = .10. 


TABLE IX. 
RESULTS FOR ROASTS A;, Bs, Cs 


t = 3,n = 10, B, = 8.2548 


A B C 


A = 1 5 34 
B 9 = 8 23 
Cc 5 2 = 33 


t = 3,n = 10, Bi = 6.6647 
xi = 0.32 (0.32) 
Sig. level of treat. = .007. 


A B C 


A - 7 6 27 
B 3 = 8 29 
C 4 2 — 34 


x} = 2.33 (2.31) 
Sig. level of treat. = .21. 


TABLE X. 
RESULTS FOR ROASTS As, Ba, Cs 
| | 
A B Cc ari 
A | | 6 sia 
4 | 32 
| | 6 27 


t = 3,n = 10, B, = 8.6191 
xi = 0.55 (0.56) 
Sig. level of treat. = .50. 
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The level of significance is approximately .6 indicating good agreement 
of the model and the observations. 

When the experiment is regrouped by sets of roasts, A, , B,, C3 --- ; 
A;,B;,(C; , the following tables and analyses are obtained. 


TABLE XI. 
RESULTS FOR ROASTS 45, Bs, Cs 


A B C 


A som 3 3 34 
B 7 saa 7 26 
C 7 3 — 30 


t = 3,n = 10, B, = 8.0688 x? = 0.51 (0.51) 
Sig. level of treat. = .13. 


For this grouping, —2 In A,; = 3.90 and, with 5 degrees of freedom, 
this value of x’ is at the .57 level of significance. 


Figure I shows the observed cumulative frequencies of values of 
xi plotted against the probability of a larger x{ as obtained using tables 


FIGURE | 
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of x’. The values of x; computed are clearly in line with what would 
be expected for a sample of ten observations on a xi-variate. These 
results are indicative of the appropriateness of the model for this ex- 
periment. We note that there is some dependence among the ten 
computed values of x; that are used in Figure I. 

(ii) Experiment 2. The Effect of Color Set on Preference, Red 
Color, and Yellow Color for Stayman Apples. G. E. Mattus, Virginia 
Agricultural Experiment Station. 

A commercial hormone preparation called ‘color set’ is used to 
delay the dropping of mature apples from the trees. The preparation 
is applied as a spray a short time before harvest to prolong the picking 
period. It was thought that the color set produced some change in the 
surface color of the apples and in particular that it increased the in- 
tensity of the yellow tones and also created the appearance of more red 
color. The following experiment was conducted to investigate these 
possibilities. 

Two groups of judges were asked to rank representative plates of 
Stayman apples for preference, red color, and yellow color in a paired 
comparisons experiment. The treatments were A (check or no color 
set), B (q p.p.m. of color set), and C (m p.p.m. of color set) with g > m. 
Seven judges, who were selected from the stenographic and laboratory 
staff, were used in one group and eight horticulturists were used in a 
second group’. The rankings of the judges in each group were pooled 
making totals of 21 and 24 repetitions of the paired comparison experi- 
ment. 

The results of tests of goodness of fit are summarized in two-way 
tables as they were for Experiment 1. Since the values of n are 21 and 
24, values of ps , Pe , and pc were obtained by direct calculation 


TABLE XII. TABLE XIII. 
COLOR PREFERENCE—GP. 1 COLOR PREFERENCE—GP. 2 
A B Cc A B C 

A - 12 6 66 A - 17 8 71 

B 9 - 3 72 B 7 - 8 81 

C 15 18 - 51 C 16 16 - 64 
t =3,n = 21, B, = 15.5184, = 3,n = 24, B, = 19.8594, 
x? = 0.43. x? =3 1.38. 
Sig. level of treat. < .005. Sig. level of treat. = .016. 


‘The results of a third group composed of statisticians have been deleted as unreliable for they 
showed extreme disagreement on color preferences. 
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(cf [4]). Here again the values of x” are not suggestive of poor agreement 
of the model with the data. One value, that in Table XVII, is rather 
large but this may be partly due to small cell frequencies. 


TABLE XIV. TABLE XV. 
MOST RED COLOR—GP. 1 MOST RED COLOR—GP. 2 
A B C 2r; A B C zr 
A - 10 8 66 A - 13 7 76 
B il - 13 60 B 11 = 12 73 
Cc 13 8 _ 63 C 17 12 -_ 67 
| 
t = 3,n = 21, B, = 18.7155, t = 3, n = 24, B, = 21.1635, 
x? = 1.30. x? = 2.11. 
Sig. level of treat. = .57. Sig. level of treat. = .31. 
TABLE XVI. TABLE XVII. 
MOST YELLOW COLOR—GP. 1 MOST YELLOW COLOR—GP. 2 
A B C 2r; A B C 2r; 
A _- 11 5 68 A _ 14 0 82 
B 10 = 5 69 B 10 - 4 82 
Cc 16 16 52 C 24 20 52 
t = 3,n = 21, B, = 16.3254, t = 3,n = 24, Bi = 13.2043, 
x? = 0.01. x? = 6.58. 
Sig. level of treat. < .005. Sig. level of treat. < .005. 


(iii) Experiment 3. The Effect of Monosodium Glutamate (M.S.G.) 
on the Flavor of Applesauce and Dehydrated Apple Slices. Lyle L. 
Davis, Virginia Agricultural Experiment Station. 

A fairly extensive experiment was conducted on the effect of various 
amounts of monosodium glutamate in applesauce and dehydrated sliced 
apples. Monosodium glutamate is frequently used as a taste stimulant. 
Much of the experiment was conducted by the method of paired com- 
parisons and evaluated by means of the method of the rank analysis of 
incomplete block designs. The score cards were immediately 
summarized by computing treatment sums of ranks and unfortunately 
the greater part of these cards were then destroyed. 

The following tables show results for only a part of the experiment. 
The treatments may be taken to be A (u p.p.m. of m.s.g.), B (v p.p.m. 
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of m.s.g.), C (w p.p.m. of m.s.g.) and D (check or no m.s.g.). For all 
sections of data presented u < v < w but for the four sections wu, v, and 
w are not identical concentrations. 

In each group each judge performed just one repetition of the paired 
comparisons design. The results of all judges in each group were 
pooled for the analysis and the tables below show the cell frequencies 
for the groups together with values of B, and x’, the latter for tests of 
goodness of fit.  x° now has three degrees of freedom. Judges recorded 
that sample of a pair that they preferred in each case. 


TABLE XVIII. TABLE XIX. 
M.S.G. IN APPLESAUCCE—GP. 1 M.S.G. IN APPLESAUCE—GP. 2 
A B | | D | 2 A} 
16 B| 4 -| 
Pie ie) 4) Di 
| 
t= 4,n = 4, B, = 5.4139, x2 = 6.94. t = 4,n = 6, B, = 8.1788, x3 = 0.76 
Sig. level of treat. = .059. Sig. level of treat. = .010. 
TABLE XX. TABLE 
M.S.G. IN APPLESAUCE—GP. 3 M.8.G. IN DEHYDRATED APPLE SLICES 
| 
A|B|c|D | Bj | D | 
A} a | -| 4 | 4 | 0 | a8 


t=4,n = 4, By = 7.115%. x2 = 0.54. t = 4,n = 7, B, = 11.2387, x3 = 6.32. 
Sig. level of treat. = .986. Sig. level of treat. = .104. 


Here two values of x; are approaching the 5% significance level 
(7.82) and two are quite small. Again, the large values appear where 
zero cell frequencies occur and in such cases x’ is likely to be a poorer 
approximation to the distribution of the appropriate likelihood ratio 
statistic than when small cell frequencies do not occur. Further difficulty 
may be introduced due to tffe small values of n. The mean value of 
x3 is expected to be 3 and we have an average value of 3.5 for this 
experiment. 
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6. Concluding remarks: 


In using a test for the appropriateness of a model as we have done, it 
is customary to use a continuity correction in computing x’. This 
tends to diminish the computed value of x° and hence to reduce its 
apparent significance. We have not used any such corrections in the 
values computed and they may be biased in such a way as to emphasize 
apparent departures from the model. However, in our data, twenty 
values of x* have been computed and only one is significant at the 5 per 
cent level of significance. 

It may be argued that the data presented are too scanty to demon- 
strate the applicability of the model for paired comparisons in general. 
The major objective of this paper is to present procedures for testing 
the models and to illustrate the applications of the procedures. As a 
secondary objective we have presented additional experimental results 
with a view to showing enough data to provide some assurance that the 
model may be appropriate for experiments involving subjective judg- 
ments. It has recently come to our attention that J. W. Hopkins has 
independently supervised extensive experiments to provide data for 
testing the appropriateness of the models for experiments on basic 
taste sensations. The results of his work have been reported [3] and 
he did not find any indication of serious departures of observations from 
the observed frequencies expected under the model. 

The author would like to express his appreciation to Lyle L. Davis, 
G. E. Mattus, C. M. Kincaid and H. R. Thomas for the use of their 
experimantal results. 
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INCOMPLETE BLOCK RANK ANALYSIS: SOME 
TASTE TEST RESULTS’ 


J. W. Hopkins 
Division of Applied Biology, National Research Council of Canada 


INTRODUCTION 


Bradley, Terry and co-workers (2, 3, 9) have put forward a solution 
of the problem of incomplete block rank analysis, and have applied it 
to the design and analysis of sensory difference tests of flavor and 
palatability. In its simplest form of paired comparisons in balanced 
incomplete blocks of 2, their method rests upon two basic assumptions. 
The first of these is that all rankings are statistically independent. The 
second is that if aliquots of m test samples are appraised by each of n 
subjects, there exist n sets of non-negative real numbers 7; , m2; , °° 
(J = 1, 2, --+ termed “ratings” and each satisfying the 
condition >>;7;; = 1, such that for all pairs of samples h and i, the 
probability of the jth subject ranking 7 above h in a single paired com- 
parison is specified by 2;;/(,; + 7;,). 

Appropriateness of this conceptual model in any specific instance 
is a question of fact. Technological interpretation of subjective test 
results also involves another factual consideration, namely the extent 
of inherent differences in the z;; characteristic of individual subjects j. 


Three series of experimental results bearing upon these points are 
reported below. 


EXPERIMENTATION 


Eleven subjects were recruited from laboratory staff. A, B, C, G 
and H were women; D, /, F, J, J and K men. All had previously 
participated in taste tests. 

In Expt. I, 6 subjects, A, B, C, D, E and F each made 20 independent 
replicate rankings of relative sweetness of the two test aliquots in each 
of 12 pairs of aqueous sugar solutions. Six of these comprised the 
possible pairings of solutions S1, S2, S3 and S4 containing respectively 
2.16, 2.32, 2.49 and 2.67% sucrose by weight in tap water; the other 6 
comprised the possible pairings of solutions L1, L2, L3 and L4 having 
6.90, 7.25, 7.62 and 8.00% lactose by weight. These concentrations 
increase geometrically by factors of 1.073 for S and 1.051 for L. Pre- 


liminary trials indicated that they provided suitable degrees of taste 
contrast. 


1N.R.C. Publication No. 3357. 
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TABLE TI. 
Individuals’ Recorded Frequency f of Specified Rankings in 20 Appraisals, and Corresponding 
Bradley-Terry expectations f’ 


Specified A B D 1 F Total 
i3.2>1 16 | 15.8 | 17 | 16.3 | 15 | 15.6 | 12 | 13.8 | 16 | 14.3 | 19 | 16.6 | 94 
3>1 17 | 17.9 | 18 | 17.8 | 19 | 18.8 | 19 | 18.0 | 13 | 16.2 |] 18 | 19.5 | 105 
4>1 20 | 18.8 | 19 | 19.8 | 20 | 19.8 | 20 | 19.3 | 18 | 18.5 | 20 | 19.9 | 118 
3>2 13 | 13.9 | 14 13.0 | 15 | 16.4 | 16 | 16.1 | 16 | 12.6 | 19 | 17.9] 92 
4>2 17 | 16.2 | 19 | 19.3 | 20 | 19.2 | 17 | 18.6 | 17 | 16.6 | 19 | 19.4 | 109 
4>3 11 | 13.0 | 20} 18.8 | 16 | 17.0 | 16 | 15.2 | 15 | 14.9] 15] 15.7] 93 
L.2>1 12] 11.0] 14] 15.0 | 15 | 15.7 | 14] 14.4 15 | 14.8] 16] 16.6] 86 
3>1 9 | 10.6 | 17] 16.6 | 19 | 17.6 | 13 | 13.7 | 16 | 16.8 | 17 | 17.3 | 92 
a> 1 16 | 15.4 | 19 | 18.4 | 18 | 18.6 | 17] 16.0 | 18 | 17.3 | 19 | 18.9 | 107 
3>2 10; 9.6] 12) 12.3] 11] 13.3] 8] 9.1] 12] 13.0] 11] 11.2] 64 
4>2 15 | 14.6 | 15 | 15.7 | 17 | 15.6 | 13 | 12.1] 15 | 13.8 | 15 |] 15.5 | 90 
4>3 14] 15.0 | 14] 13.9 | 12] 12.8] 11] 12.9] 9] 10.9] 15] 14.6] 75 
Specified B E G H I J Total 
ranking — f 
ILS.2 > 1 It | 12.7] 15] 15.8] 9] 11.9] 9] 10.3] 13 | 11.8 | 12] 12.1 69 
3>1 18 | 16.8 | 17 | 16.9 | 17 | 15.4] 16] 14.6] 15 | 15.5 | 17] 15.2 | 100 
4>1 20 | 19.7 | 20 | 19.4 | 20 | 18.7 | 16 | 16.1 | 17 | 17.6 | 15] 16.8 | 108 
3>2 14 | 15.0 | 12 | 12.0 | 12 | 13.8 | 14) 14.4] 16] 14.1] 14] 13.5] 82 
4>2 19 | 19.6 | 17 | 17.8 | 17 | 18.1 | 15 |] 15.9 |] 16 | 16.7] 15] 15.6 | 99 
4>3 19 | 18.7 | 17 | 16.9 | 16 | 16.2 | 13 | 12.0] 15 | 13.7] 15 | 12.6] 95 
SM.2>1 13 | 15.7] 10} 11.4] 9] 11.2 | 14] 15.0 | 18 | 16.8 | 13 | 12.2 | 77 
3>1 19 | 16.9 | 17 | 14.2 | 12] 10.5 | 19 | 18.0 | 19 | 18.6 | 12 | 13.2] 98 
4>1 20 | 19.6 | 16 | 16.4 | 12 | 11.2 | 19] 18.9 | 18 | 19.5 | 17 | 16.6 | 102 
3>2 10 12.0] 11] 12.9] 8] 9.2] 13] 15.0] 15] 14.3] 12] 11.0] 69 
4>2 18 | 18.5 |] 15 |] 15.5 9} 10.0 18 | 17.0 | 18 | 17.6 | 15] 15.2 93 
4>3 18 | 17.9 | 14 | 13.1 | 11 | 10.8 | 12} 13.1] 16 | 14.9 | 14] 14.3] 85 
| 
Specified B D E G H kK Total 
f ig f f | id f f f’ f 
Il.W>X 4 5.4] 15] 13.2] 5 5.0) 8| 9.5] 8] 8.2] 10]10.3] 50 
16 | 14.2 | 17 | 18.2 3 3.1 6.5] 7] 5.56] 5] 6.0] 57 
W>Z 8.3] 12|12.6] 9] 88] 7] 8.0] 9] 10.3] 10] 8.7] 55 
4 17 | 17.3] 16.8] 8] 7.1 8} 7.0} 6] 7.1 5| 5.8] 61 
X>Z 12 | 13.1 | 1) 9.4) 13 | 14.0 6 8.5) 13 | 12.1 9 8.5 64 
Y>Z 6 4.5 2 2.9 | 17 | 16.2 | 15 | 11.5] 15 | 14.7] 11 |] 12.8 66 


In Expt. II, subjects B, E, G, H, I and J also each made 20 in- 
dependent replicate rankings of relative sweetness of the two test 
aliquots in each of 12 pairs. Six of these were identical with the fore- 
going pairs of S1, S2, S3 and S4; the other 6 were pairings of solutions 
SM1, SM2, SM3 and SM4 derived from them by addition to each of 
0.003% quinine sulfate, 0.15% sodium chloride and 0.21% tartaric 
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acid to provide uniform accompanying bitter-salt-sour taste stimuli. 

For Expt. IIT, four qualitatively different test solutions, W, X, Y 
and Z were used. These were arrived at from a basic mixture of 0.625% 
sucrose, 0.075% sodium chloride, 0.02675% tartaric acid and 0.00015% 
quinine sulfate by weight, again in tap water, by reducing the con- 
centration of each of these constituents in turn by 1/4. Subjects 
B, D, E, G, H and K then made 20 independent preference rankings of 
these solutions in all 6 pairings. 

In all experiments, each subject ranked six pairs of aliquots at 
daily mid-morning and mid-afternoon tasting sessions. The order of 
presentation of test pairs and of aliquots within pairs was objectively 
randomized subject to 2 restrictions. (i) Each pair of aliquots was 
presented 10 times in each of the sequences Ai and th. (ii) In Expts. 
[ and II, daily morning and afternoon sessions 1 and 2, 3 and 4 etc. 
together comprised all 12 test pairs, while single sessions each comprised 
three S and three L or three S and three SM. All aliquots were brought 
to 70°F. for tasting in a laboratory conditioned to an air temperature 
of 72°F. and a relative humidity of 50%. 


ANALYSIS OF RESULTS 
Repeatability 


Table I lists the frequency with which individuals recorded specified 
rankings in their 20 appraisals of each aliquot pair, S2 > 1 denoting 
ranking of S2 above S1, ete. Fig. 1 portrays the repetitive fluctuations 
in the 6-subject 6-comparison aggregate frequencies of these rankings 
recorded for each series of solutions. (In Expt. III each session, and 
in Expts. I and II the morning and afternoon sessions of each day 
together constituted a complete replicate). No evident progressive 
trend in these totals is visible. A numerical trend test was made by 
ordering the 20 replicate total frequencies of all 12, 12 and 6 specified 
rankings in each experiment, assigning the appropriate Fisher-Yates 
(7) normalized rank score to these ordered frequencies, and computing 
the variance associated and not associated with any linear trend in the 
resulting scores. (Under conditions of statistical stability these scores 
should of course have the characteristics of random normal deviates). 
Variance ratios of F(1, 18) = 1.56, 1.44 and 1.13 for Expts. I, II and 
III respectively and of I'(3, 54) = 1.45 in the aggregate resulted. As 
the 5% points of F(1, 18) and of (3,54) are 4.41 and 2.78, this analysis 
is likewise indicative of no progressive group trend. 

As an index of the statistical stability of individuals’ appraisals, 
Cochran’s (5) Q for matched samples was computed from the replicate 
frequencies of specified rankings by each subject in each experiment. 
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The resulting Q are brought together in Table II. In the absence of 
heterogeneity, these individual Q should be distributed approximately 
as x’ with 19 df. (5% point 30.1), their sub-total for each experiment 
as x’ with 114 d.f. (5% point 140 approx.) and their total for all 3 
experiments as x” with 342 d.f. (5% point 386 approx.). In fact only 
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FIGURE 1. REPETITIVE FREQUENCY OF SPECIFIED RANKINGS IN EXPERIMENTS 

I, 1 AND III. S, SUCROSE ALONE; L, LACTOSE ALONE; SM, SUCROSE WITH OTHER 

PRIMARY TASTE STIMULI; W, X, Y, Z, QUALITATIVELY DIFFERENT TASTE 
MIXTURES. 
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TABLE II. 
Indices Q of Discrepancies Between Replicates in Recorded Frequency of Specified Rankings by 
Individual Subjects 


Experiment I Iexperiment II Experiment III 
Subject Q (19 d.f.) Subject Q (19 d.f.) Subject Q (19 df.) 
A 29.7 B 17.6 B 15.5 
B 28.7 E 23.3 D 18.8 
C 22.8 G 19.6 E 14.6 
D 23.2 H 9/2 G 18.7 
E 18.4 I 26.1 H 18.3 
F 23.8 J 33.2 K 18.7 

Total 146.6 Total 129.0 Total 104.6 


one of the 18 individual Q listed, that for J in Expt. II, exceeds its 
5% point. Subtotal Q for Expt. I attains approximately its 2% point, 
but as this is the largest of 3 subtotals it is suggestive rather than 
clearly indicative of correlated fluctuations in the rankings recorded in 
different replicates by the same subject. The other sub-total Q are 
entirely consistent with statistical independence of all rankings in 
Expts. II and Iil. 


Bradley-Terry Rank Analysis 


In Bradley and Terry’s analysis (2), maximum likelihood estimates 
pi; of the “ratings” 2,; from s replicate rankings of all pairings of m 
test samples are specified by the equations: 


+ pi) = 0 = 1,2, m) 
Pii hei 
=1 
in which 
k=s 


Tisjk being the rank assigned by the jth subject to the 7th test sample 
in the kth replicate paired comparison of h and 7. Approximate solution 
by iteration of these equations for the rank frequencies recorded in 
Table I gave the p,, listed in Table III, and hence the Bradley-Terry 
expectations f{; = 20p,;/(pa; + p;;) also listed in Table I. 
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TABLE III. 

Maximum Likelihood Estimates of Bradley-Terry ‘‘Rating’’ Parameters for Individual Subjects 
Test A B C D 
item 

$5.1 034 007 .009 .024 050 005 
2 129 031 .032 .053 124 024 
3 294 058 .143 .221 212 208 
4 543 904 .816 . 702 614 763 
150 050 .040 . 103 065 034 
2 184 151 . 148 . 266 183 168 
3 168 242 .291 . 223 343 216 
4 498 557 .521 .408 409 582 
Test B E G H I J 
item 
ILS. 1 012 025 .050 113 | .075 091 
2 021 093 .074 120 109 138 
3 063 138 . 166 .305 259 285 
4 904 744 .710 .462 557 486 
IL.SM. 1 018 107 .214 .033 017 106 
2 065 142 . 275 .100 090 167 
3 098 260 . 236 .300 228 206 
4 819 491 .275 .567 665 521 
Test B D Dy G H K 
item 

Ill. W #82 454 .094 176 167 180 

X 488 236 .278 196 240 171 
= 075 045 . 509 362 435 416 
Z 255 265 119 266 158 233 


Goodness of fit of the observed and expected frequencies was then 
tested as suggested by Bradley (4) by computing for each subject and 
test series the index of discrepancy }-(f — f’)?/f', this summation 
comprising in each instance the complementary frequencies of rankings 
« >handh>vz. The results are shown in Table IV. On the Bradley- 
Terry assumptions, individual indices should be distributed approxi- 
mately as x” for 3 d.f. (65% point 7.82), their 6-subject subtotals approxi- 
mately as x” for 18 d.f. (5% point 28.9) and their total for all three 
experiments, here 84.2, approximately as x’ with 90 d.f. Clearly the 
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TABLE IV. 
x? (3 df.) from Individual Tests of Goodness of Fit of Recorded and Expected Rank Frequencies 


Experiment Experiment I-xperiment III 
Subject | I, Subject | | SM Subject | WXYZ 
A 299 | 1.11 | B | 2.72 | 5.17 B 2.24 
B 5.04 | 0.80 | E 1.14 | 3.40! D 2.53 
2.21 | 323 G | 534 | 20 | E 0.63 
E 3.82 |1.7| 1 199 | 5.78 | 4H 1.36 
F 6.83 | 0.46 | J | 3.49 | 075 | K | 1.58 
! | | 
Total | 25.06 | 9.01 | Total | 15.99 | 19.72 | Total | 14.38 
(18 (8 df.) | df) | 
| | | 


numerical values listed are individually and collectively indicative of 
excellent agreement between the observed and expected rank frequencies, 
notwithstanding the fact that several of the expectations f’ complemen- 
tary to those given in Table I are less than unity. 


Difference Between Subjects 


Homogeneity of individual reaction to the same paired comparison 
was tested directly by arraying the frequencies of ranking 7 > h and 
h > 7 recorded for each such comparison by the 6 participating subjects 
in the form of a2 X 6 contingency table, computing x’ from these tables, 
and adding the results for each test series 7.S, J.L etc. Table V shows 
these aggregate x”. Here again some of the expectations were small. 
In this instance however Haldane’s (8) formulae permitted calculation 
of the exact moments of the conditional distribution of x’ for the 
marginal totals of each 2 X 6 table. The resulting aggregate K,(x’) 
and K,(x’) are accordingly also listed. 

For the series J.S, 7.L and J/.S, (90 is 115.1, exceeding 
its expectation of 90.9 by 24.2, which is 1.88 times the corresponding 
standard deviation of (166.88)'’*, namely 12.9. On the basis of upper 
tail probabilities this is suggestive, but only suggestive, of some real 
differences in individual discrimination of intensities of sweetness. 
There are evident differences between some of the corresponding total 
rank frequencies for 7.S and //.S listed in the last column of Table I, 
but as the two tests were made at different times experimental factors 
(e.g. differences in tap water) may have contributed to these 


crepancies. 
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Individual differences in respect of the 77.SM comparisons were 
much more pronounced, and gave an aggregate x” (30 d.f.) exceeding 
its expectation by 6.3 times the corresponding s.d. Indications from 
Tables I and III are that the discriminatory acuity for sweetness of 
subjects E and G was notably reduced by presence of the other three 
tastes; that of B and J was seemingly little affected, while that of H 
and J was if anything increased. 

Individual differences in respect of comparisons III were still more 
pronounced, the aggregate x’ (30 d.f.) of 124.5 for these shown in 
Table V exceeding its expectation by 12.3 s.d. No two of the six 


TABLE V. 
x? from Tests of Homogeneity of Frequency of Specified Rankings by Individual Subjects 
Test x? (30 d.f.) Ki (x?) Kz (x?) 
items 
I. 8. 45.1 30.3 52.18 
35.1 30.3 57.58 
II. S. 34.9 30.3 57.12 
II. SM. 78.2 30.3 57.67 
III. WXYZ. 124.5 30.3 58.44 


participants rated all 4 test flavors in the same order. Y was rated 
first by 4 subjects but last by the other 2. X and W were both rated 
last by 2 subjects but first by another. 


DISCUSSION 


Repeatability and mutual independence of psychophysical reactions 
can seldom be taken for granted. In the present experiments, however, 
these conditions seem to have been reasonably approximated. 

The Bradley-Terry model implies (3) that the probability of 7 
being ranked above h in a single paired comparison by the jth subject 
may be specified by the logistic function 
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in which y = In (x,;/7,;) may be regarded as a measure of the degree 
of sensation evoked by the taste contrast (i, h). The generally good 
agreement of observed and hypothesized rank frequencies might then 
be regarded as consistent with existence of a specifically ordered con- 
tinuum of degrees of sensation characteristic of each individual. For 
moderate frequencies the sigmoidal graph of a logistic is virtually in- 
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distinguishable from that of an integrated Gaussian function. Hence 
when this model is appropriate and stimulus differences are expressible 
in terms of a measured composition or process variable, knowledge of 
the appropriate transform f(z, h) = y would permit repeated paired 
comparisons of a series of qualitatively similar test sample s to be 
regarded as organoleptic assays treatable statistically by the standard 
methods of logit or probit analysis (6). 

In the foregoing experiments individual subjects reacted differentially 
to all but possibly the simplest flavor contrasts. Bradley and Terry 
(2) have noted that in these circumstances the ‘‘ratings”’ 7, ; are functions 
of the subject as well as of the test material. Technological applications 
may then require numerical prediction of population reactions from those 
of small groups of laboratory subjects. 

Opinion regarding commensurability of “degrees of preference’’ on 
a subjective continuum is currently divided (1, 10). It is therefore 
perhaps noteworthy that in these trials individuals’ preference rankings 
of the qualitatively different test solutions were as consistent with the 
Bradley-Terry conceptual model as were the rankings of single flavor 
intensities. 
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A NOTE ON MISSING-PLOT VALUES 


J. A. NELDER 
National Vegetable Research Station, Wellesbourne, Warwick, England 


This note owes its origin to the discussion in Biometrics between 
Professor Snedecor and Dr. Williams on Query 96 (Biometrics Vol. 8 
No. 4). In‘his original answer Professor Snedecor repeats what has 
often been written when he says ‘The value of X is not intended as an 
estimate of the missing datum’. Nevertheless whether or not X is 
intended to be an estimate of the missing datum, it zs an estimate of 
the missing datum, and an unbiased one where the mathematical model 
used is true. Thus, in the randomised block design, a commonly used 
model is that in which x;; = m + 8; + 7; + €;; where z;; is the observa- 
tion in the jth treatment in the 7th block, m is the mean, 8 and 7 are 
(fixed) block and treatment constants respectively (subject to the usual 
constraints >> 8 = 0, >> r = 0) and e«,; are the random residuals, 
mutually uncorrelated. If it is further assumed that the variance of 
€;; is constant, i.e. independent of 7 and j, then the formula for the 
missing plot value (y) in a randomised block experiment is y = 
(rB + tT — G)/[(r — 1)(t — 1)], where r = no. of replicates, ¢ = no. of 
treatments, B is the total of the remaining values in the block with the 
missing plot, 7’ the similar total for the treatment, and G is the total 
of all existing values. 

This value is the one obtained in a full least-squares analysis of the 
complete experiment, and its expectation is m + 8; + 7; , i.e. the mean 
value of the missing datum. Its variance is a minimum among linear 
unbiased estimates by the familiar properties of least squares. The 
variance of y for the randomised block design is given by 

‘, 
var (y) = 


and for the r X r Latin square by 


ew = - 2° 


It is important to note that if the conditions postulated above do not 
hold, for example, if the variance of the e;; is a function of 7 and j, then 
the missing plot value given above is no longer correct, and does not lead 
to the correct analysis. 
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An ‘impossible’ value of y (such as a negative yield) or an un- 
reasonably discrepant one can arise through sampling error, or it can 
arise through a failure of the mathematical model, and it may be of 
some importance to distinguish which is the cause in any particular 
instance. On the assumption that the missing plot has occurred at 
random in the experiment, we could use the estimated variance of the 
missing plot value to form a confidence interval in the usual way. If 
such a confidence interval for an essentially positive quantity contains 
no positive values, this is evidence against the negative value’s being 
due to sampling error, and so in favour of its being due to a failure of 
the model. The practical value of this test is limited, however, by our 
knowing little about the expected value of the missing-plot value, other 
than that it must be positive if it is a yield or a count, etc. In Query 96, 
the test gives —6.64 + 8.10 as the 95 per cent. confidence interval for 
y. Hence positive real values for y are not ruled out, and yet it is 
fairly certain from other evidence, as Dr. Williams points out, that a 
multiplicative model is much nearer the truth here than an additive 
one. Even when the test fails to discredit the sampling error hypo- 
thesis it is still hardly sufficient to suppose that the estimation of the 
error variance and the tests of significance will be necessarily valid. 
In Query 96, the variance of the catches in the treatment with the 
highest mean (B) is over 8 times the variance of the catches in the 
treatment with the lowest mean (A). The use of a common error 
variance for all comparisons will give considerably biased standard 
errors for some of the comparisons. Here, it happens not to matter, 
but this is nothing to take comfort from. 

The analysis of variance receives rough treatment from some of its 
practitioners, and it is one of its merits that it can stand much of this 
without serious breakdown. However, it is as well not to try it too 
hard. An ‘impossible’ missing plot value is a danger signal, and should 
be taken heed of. 


4. 


A NOTE ON TRUNCATED POISSON DISTRIBUTIONS 


P. G. Moore 
University College, London and Princeton University 


1. In an observed Pcisson distribution we have counts of the number 
of times that a certain event has occurred once, twice and so on in a set 
of N readings. From theoretical considerations we obtain the probability 
of just r occurrences of the event to be 


+r =0,1,2,--- 


where d is the mean of the distribution and thus represents the average 
number of events to be expected per reading. An observed set of data, 
although Poisson in form, may exhibit various types of truncation due 
to the whole of the data not being available. If n, be the observed 
frequency of just r occurrences then the main types of truncation are 


(1) n, for r > s not observed. 

(2) n, for r < k not observed. 

(3) combination of (1) and (2). 

(4) n, for r = v tor = w inclusive not observed. 


Truncations of Type 1 arise when there is difficulty in counting high 
numbers owing to inability to distinguish each individual or when the 
counting apparatus gives trouble for high counts. Type 2 occurs when 
it is desired to estimate the mean frequency of some event and only 
cases of more than k occurrences are reported. Type 3 occurs sometimes 
in botanical work where only quadrats with at least one individual in 
them are retained and high densities are difficult to count. Type 4 
could occur when the observations have been muddled in some way. 

Type 1 has been considered by Tippett [11] who derived the maximum 
likelihood solution, by Bliss [1] who gave tables to facilitate the solution 
of the equation and by Moore [6] who suggested a simple form of esti- 
mate needing no tables but having a larger standard error. Type 2, in 
the case where k is equal to zero, has been studied by David and Johnson 
[2] who gave the maximum likelihood solution, by Plackett [2] who 
gives a solution similar to that of Moore for Type 1 and Rider [9] who 
uses the first two incomplete moments. Cohen [3] has provided maxi- 
mum likelihood estimates similar to those of Tippett for the cases of 
Types 2 — 4 and indicated how to solve the resulting equations by an 
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iterative or interpolatory procedure using a table of the Poisson dis- 
tribution. In this note we will show that the type of estimate given in 
our earlier paper may be applied to any of Types 1 to 4 truncations by 
merely making a suitable change of the limits of summation. This 
quick form of estimate might be useful in any of the following situations 


(a) when no extensive table of the Poisson distribution is available, 

(b) when it is desired to obtain an estimate quickly without much 
computation, 

(c) as a first approximation to use in the full maximum likelihood 
solutions of Tippett and Cohen. 


It is, of course, only to be expected that the standard error of such an 
estimate will be larger than for the corresponding maximum likelihood 
solution. 


2. For a Poisson distribution, truncated as in Type 1, we may write 
the identity 


and this suggests using as an estimate of \ the expression 
Type 1 1 = Dm,/ dn, 
r=0 r=0 


A similar form of argument for the other three types of truncation will 
give us as our respective estimators 


Type2 }, = rn, n, 
r=k+2 r=k+1 
Type3 }; = / n, 
r=k+2 r=k+1 
9-2 o 
Type 4 {Som + + =. 
r=0 r=wt+2 r=0 r=wt+i 


Thus the only quantities used in the various estimators are r and rn, . 


3. As an example we will consider the following data due to Scrase 
{10}. 


oo 


Total 


r 0 1 2 3 4 | | as 
n, | 23 | 56 | 88 | 95 | 73 | 40 | 17 | 5 | 3 | 400 ae 
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The problem is connected with the number of dust nuclei in the air and 
the data give the frequency distribution of the number of drops in a 
small volume of air that fall on to a stage in a chamber containing 
moisture and filtered air. The estimates obtained using the various 
formulae are tabulated in Table 1 where the estimate is given after the 


TABLE 1. 
Type 1 s=5 3.088 s=6 3.012 s=7 2.963 
Type 2 k=0 2.955 k=1 2.922 k =2 2.803 
Type 3 k=0 k=0 k=1 
2.997 3.054 2.970 
8=7 s=6 s=7 
Type 4 v=4 
3.000 
w=5 


differing truncation details. For the full distribution the mean is 2.925 
but Scrase is of the opinion that this mean is slightly high in that a 
number of zero counts were wrongly rejected as being due to the 
apparatus not working. 

4. There are two ways of looking at the standard errors of the 
estimates. To show this let NV be the actual number of observations 
counted, that is N less those individuals that have been truncated. 
Then in obtaining the standard error we can regard either N or alter- 
natively N as being constant from sample to sample. The first situation 
was used by Tippett and Moore whilst the second, which appears to 
lead to simpler results, was used by David and Plackett. It is clear 
that to lay down hard and fast rules as to which is the correct one to 
use is impossible and situations of each type arise. For instance in 
particle counting where cells with high frequencies are not counted or 
in accident returns for some fixed population submitted to an authority, 
the situation with N constant seems plausible. On the other hand if 
the population is varying and we have just the data concerning actual 
accidents we would regard N rather than N as fixed. 

Assuming N to be fixed we can obtain asymptotic expressions for 
the standard errors of truncations of Types | and 2 fairly easily by 
the use of expectation techniques. ‘They simplify to 
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r=0 r=0 


o 2 2 
r=k+2 r=k+2 r=k+1 
The variance of \, has been discussed by Moore, but for \, we append a 
few values in Table 2 when k is equal to zero, one and two. To compare 
the case of k equal to zero with the figures given by Plackett we should 
multiply by (1 — ,)/N since his total is for the non-truncated portion 
only. On doing this we find our results comparable with his. In 
situations where the whole distribution is available and X is estimated 
from the mean, N Var (mean) is simply \. This enables us to see from 
Table 2 the loss of accuracy due to the truncation of the data. 


TABLE 2. 
N Var 
k=0 k=1 k=2 
1.0 3.0846 11.8376 52.4339 
1.5 3.4089 10.0824 32.8254 
2.0 3.6256 9.1150 24.8133 
2.5 3.8914 8.5522 20.6345 
3.0 4.1250 8.0800 17.8645 
3.5 4.3837 7.7107 15.8528 
4.0 4.6772 7.4367 14.2889 
5.0 5.3713 7.1754 12.0162 
6.0 6.1941 7.2961 10.6001 


5. The foregoing procedure may, with some small modification, be 
used for the case of truncation of the binomial population 


(q + p)” 


where n is known but p is unknown and the number of zeros is missing. 
Cases of this kind arise, for example, in genetical work where it is 
desired to estimate the proportion of sibs who inherit some abnormality 
if the parents possess a certain gene. The presence of the gene is, how- 
ever, only known if at least one of the sibs possesses the character in 
question. Hence the number of zeros in sibships of a given size n 
cannot be known and we must use the remaining information by itself 
to estimate p. Fisher [5] gave the maximum likelihood solution and 
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Finney [4] gave tables to reduce the solution to a fairly rapid iterative 
procedure. Following our previous method we are led to using 


n rn, n-1 


as an estimator of p/q and from this equation we obtain 


n rm, n-1 n rn, \ 
as an estimator for p. 
To illustrate this procedure we consider the following data, due to 
Pearson [7], which give the number of sibs with human albinism in 


families of size five. Substituting in the formula given above we get 
0.326 as our estimate of the unknown probability. 


r 1 2 3 4 5 
Nr 25 23 10 1 1 
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QUERIES 
GrorGe W. SNEDEcoR, Editor 


QUERY: Data on the calcium, phosphorous, and magnesium 
110 content of turnip leaves were collected in the following manner 
as reported in Southern Cooperative Series Bulletin No. 10, 
1951, “Studies of Sampling Techniques and Chemical Analyses of 
Vegetables”. ‘Duplicate [microchemical] analyses were made on each 
of four randomly-selected leaves from each of four turnip plants picked 
at random --- . Duplicate determinations were made on each ash 
solution from a particular leaf --- . The analyses of the two sets of 
ash solutions were made at different times.” 
To simplify the form of the questions to be asked, the analyses of 
variance for calcium and phosphorous are given. 


ANALAYSES OF VARIANCE 


Calcium Phosphorous 
Variance Source 

Degrees of Mean Degrees of Mean 

Freedom Square Freedom Square 
Total 60 61 
Plants 3 6.202154** 3 0.056375** 
Ashings 1 0.02945N.S. 1 0.000467N.S. 
Plants X Ashings 3 0.033569** 3 0.000664N.S. 
Leaves in Plants 12 0.605917** 12 0.035786** 
Leaves X Ashings ll 0.013968** 11 0.000935N.S. 
Duplicates 30 0.003560 31 0.000457 


**Significant at the 1 per cent level of significance. 
N.S. Not significant at either the 1 per cent level or the 5 per cent level of significance. 


It is evident that the mean square for duplicate determinations must 
have been used as an estimate of error for testing the significance of 
Plants X Ashings, otherwise this term would not have been significant 
in the calcium analysis. Would not a more appropriate test have been 
made using Leaves [in plants] x Ashings as an estimate of error for 
Plants X Ashings? 

In the case of testing for significance of Plants it is not quite so 
evident what source of variation was used as an estimate of error, but 
from the Phosphorous analysis it may be seen that of two possible 
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“contending” terms for error, Leaves in Plants and Plants X Ashings, 
each of which seems to possess the characteristics of a proper error 
term, Leaves in Plants was not used, otherwise Plants would not have 
been significant. How was it decided that Leaves in Plants would not 
be used as error for Plants? Or to make one generalized question out 
of the two, how does one select an appropriate error term in such an 
analysis as this? 


This is a troublesome problem which arises when a factor 
ANSWER: under investigation is tried over two or more other 

factors which may be regarded as random. A number of 
articles shedding light on various phases of this problem have already 
appeared in this Journal: Wilm, H. G., V:1; Satterthwaite, F. E., 
V:2; Crump, 8. L., V:2; Eisenhart, C., V:3; Cochran, W. G., V:7; and 
Ostle, Bernard, V:8 (Query). 

The choice of a proper error term in any F test of significance 
depends upon the components of variance entering the two mean 
squares. The composition of any mean square can be determined only 
after the assumptions are made as to whether the factors under study 
are random or fixed. 

In this particular case it would seem only reasonable, since the 
investigators were concerned primarily with sampling technique, and 
since they went to the trouble of specifying that both plants and leaves 
on plants were randomly selected, to assume (helped by the phraseology 
“Duplicate analyses ---”) that the two ashings are also to be regarded 
as a random sample of a large number of ashings that might result 
from repeating the experiment many times on the same leaves by taking 
at successively later dates a new random sample from each leaf to be 
ashed and then analyzed in duplicate. 

Under the assumption that all of the variables are random we may 
write the expectations of the mean squares for the various sources of 
variation as in Table 1. : 

For a case quite similar to this, but in which two of the factors were 
regarded as fixed, see Query 95, this Journal, Sept., 1952. 

The proper mean square to use as an estimate of error for testing the 
significance of any source of variation is that mean square whose 
expectation is the same as the expectation of the mean square being 
tested except for that one component which is directly due to the 
source of variation being tested. This is really a test of the null hypo- 
thesis that the component due directly to the source being tested equals 
zero, for if we set this component equal to zero then both mean squares 
have the same expectation. 
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TABLE 1. 


Expected Mean Squares of the Several Sources of Variation 


Source of Variation E(MS) 


Plants—(P) op + + + + 
Ashings—(A) > + + Sop, + 320% 
Plants X Ashings—(P4A) Op + + Sora 
Leaves in Plants (LP) op + + 
Leaves in Plants X Ashings o> + 2o7p4 
(LPA) 
Duplicates (D) oD 


To find, then, an appropriate error term for testing any particular 
mean square one looks for a mean square having the same expectation 
as has the mean square to be tested when the component due directly 
to the source being tested is set equal to zero. 

For example, the expectation of Leaves in Plants X Ashings is 
op + 2oipa. Setting op , (the component due directly to Leaves in 
Plants X Ashings) equal to zero, there remains co; . The mean square 
having expectation 5 is Duplicates. Thus the mean square for dupli- 


cates is the appropriate error term for testing differences among Leaves 
in Plants X Ashings: 


MS» 


F = 


As you suggest, the appropriate error term for testing Plants x 
Ashings mean square is Leaves in Plants X Ashings mean square. The 
expectation of Plants X Ashings is o” + 2o;p4 + 80p4 , and if the 
component due directly to Plants X Ashings (o>,) be set equal to zero, 


there remains 05 + 207», which is the expectation of Leaves in Plants X 
Ashings. For calcium: 


— MSra _ 0.033569 _ 
MSir, 0.013968 


On comparison with the tabled values / = 2.40 (with 3 and 11 
degrees of freedom) is found to have a probability of occurrence greater 
than 0.05 due to chance alone, even when no additional component 
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exists in the numerator. Plants X Ashings is therefore not to be regarded 
as significant, contrary to the analysis quoted. 

The test of significance for Plants is more complex. The expectation 
of the mean square for plants is 05 + + Sop, + + . 
Setting the component directly due to Plants («7) equal to zero the 
expectation is 05 + 2o0tr, + 8074 + 4c0zp . Examination of the 
expectations of the several mean squares of the analysis reveals that 
none of them has this expectation. However, a mean square with this 
expectation may be constructed as follows: 


MSs = + MSip — 


If we set down the expectations of the various mean squares of the 
constructed mean square and perform the specified operations we have: 


= (05 + + 8074) + (05 + + — (05 + 
= oD + + + 4oip 


Thus the constructed mean square has exactly the expectation 
necessary to serve as error for testing significances of differences among 
Plants. It is used in the same manner as any other mean square: 


_ MS>p 


F’ = 


The degrees of freedom for the constructed mean square are not 
known but may be approximated by a method due to Satterthwaite: 


(MSs)* 
(MSps)* 4 Str)” 4 Stra)” 
d.f.rs d.f.rp d.f.rps 


It has been pointed out however that the presence of the minus sign 
in the construction of MS; , especially if the degrees of freedom for the 
various mean squares are small so that the estimates might not be the 
most reliable, could possibly lead to a negative estimate of MS; in 
which case the approximating F distribution would not be too good. 
To overcome this drawback it has been proposed (Cochran) that the 
test be preformed by the calculation of: 


_ MS’ MSp + 
MS’, MSp, + MSzp’ 


= 


Fr" 


which gives two mean squares whose expectations are the same under 
the null hypothesis that o> = 0. Degrees of freedom must now be 
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approximated for both numerator and denominator of the F ratio. 


= 
(M 4 (MSzpa) 
d.f.p d.f.rpa 


Solving for both F’ and F”’ using the analysis of phosphorous content 
we find: 


~ 0.000664 + 0.035786 — 0.000935 0.035515 
(0.035515)? 

 (0.000664)* (0.035786)" (0.000935) 
3 12 ll 
0.001,261,315 _ 
0.000,106,946 


Assuming now that F’ with 3 and 12 degrees of freedom is dis- 
tributed as F we look in the Table of F and see that values as large as 
1.59 (3, 12 d.f.) occur in more than 10% of trials due to chance alone, 
so that plants are not to be regarded as differing significantly in phos- 
phorus content. 


(2) = 0.056375 + 0.000935 — 0.057310 _ 


= 0.000664 + 0.035786 ~ 0.036450 
(0.057310)? 0.003,284.436 _ 
 (0.000035)" ~ 0.001,059,459 1 = 
3 
(0.000664)? 4. (0,035786)" 0.000,106,867 
3 12 


We now have F”’ = 1.57 (3, 12 d.f.) from which we reach the same 
conclusion as before, these data are insufficient to say that differences 
in phosphorous content exist among plants, contrary again to the 
quoted analysis. 


E. F. Scuuttz, Jr. 
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QUERY: Comment on Queries Nos. 96 and 103. 
111 I find myself in partial agreement and partial disagreement with 

both querist and answerer in Query 103. The issues raised by this 
query seem important enough to warrant further discussion. The 
main point at issue is the usefulness of the logarithmic number of 
insects. This is by no means belittled by the fact that both counted 
and logarithmic values yield significance. The question is not the 
existence of insight into the data but the extent of the insight. 

We may inquire into this by examining the values in more detail. 
When we do this, the logarithmic values produce curious suggestions 
concerning the 12 trap locations—suggestions which were covered up 
in the counted values. If these suggestions correspond to real features 
in the placement of the 12 platforms, there is no doubt that the log- 
arithmic values have produced more insight. (And even if the suggested 
features are not known to be real, it is highly likely that they represent 
extra insight.) 

These suggestions may be obtained by comparing individual traps 
with the means of other traps on the same platform. (They were 
actually observed by calculating apparent interactions of the form 


(individual value)—(row mean)—(column mean) + (grand mean) 


but we shall forego these in the interests of space.) 


Reps A compared with mean of B, C, and D 
counts logarithms 
L223, 4 —47, —32, —32, —13 —.82, —.85, —.68, —.71 
5,7, 9, 11 —35, —19, (—22), —10 —.45, —.43, (—.45), —.37 
6, 8, 10, 12 —10, —13, —5, —21 —.18, —.18, —.09, —.26 


The clear distinction between the three sets of 4 replications shown in 
the differences of logarithms is dimly, if at all, perceptible in the differ- 
ences of counts. Here the logarithm has’ earned its salt—not just 
because of equal or unequal weighting—but through clarification of the 
relations of individual values. 

In passing, let me record my support of both the statement by 
Snedecor that the supplied value ‘is merely a number to be put in the 
empty space’ for convenience and the implication by Williams that the 
appearance of unusual supplied values (negative, very large, or what 
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have you) should always be taken as a warning and a basis for careful 
consideration of how the analysis used relates to reality. 

Incidentally, application of the computations of “One degree of 
freedom for nonadditivity” (Biometrics 5 (1949) 232-242)) to both the 
counts and the logarithms yields the following results: 


MS for MS for 

Counts Logarithms 
Nonadditivity 1 668 
Balance | 149 .036 


(*1 supplied value) 


The test for nonadditivity of counts is significant at 5%, while the test 
for nonadditivity of logarithms provides no slightest ground for suspicion. 
(The ratios of row MS and column MS to balance MS were each in- 
creased by a factor of about 1.5 by changing from counts to logarithms.) 
Thus an available test could have been applied to this rather restricted 
body of data to indicate the need of a transformation. 


JOHN TUKEY 


Editorial Comment: 


On receipt of this penetrating observation, I asked querist for any 
information he might have that would serve to explain the peculiar 
groupings in the data. None was forthcoming, but three more batches 
of data’ from the same traps were sent. These were examined in the 
manner suggested by Dr. Tukey. In all four batches, suspension . 
platforms 1 and 4 showed above average differences between the catch 
in trap A and the mean catch in traps B, C and D. On suspension 8, 
all differences were notably small. 

Querist sent a map of the locations of the traps but these seemed to 
have no relation to the peculiarities observed in the data. 

I am inclined to think that Dr. Tukey, by shrewd examination of the 
data, has discovered faults of construction in the traps. 
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THE BIOMETRIC SOCIETY 


Finances for 1953. On recommendation of the Finance Committee, 
the financial reports of the Society will be printed in the future in 
BIOMETRICS instead of being sent separately to each member as in 
the past. The following are the audited statements for 1953, somewhat 
condensed. 


BIOMETRICS 
Income 
From 1952, reserve from sale of No. 1 in Vols.3and7 §$ 1100.00 
Balance from 1952 books . . . ie ee 3758.21 $ 4858.21 
Member subscriptions 
| | 3992.75 5439.75 
Non-member subscriptions, 559... 3914.75 
Sale of back issues 
Vol. 3, No. 1, 121; Vol. 7, No.1,20 ...... 221.50 2593.18 
Expenses 
Overpayments, cancellations $ 77.00 
Envelopes, stationery, supplies ...........2....-. 432.90 
Institute of Statistics, Editorial Management .......... 1000 .00 
ASA, } net profit on sale of vols. 1-5. 568 .37 
Wm. Byrd Press 
Production of Biometrics, Vol.9,1953  ...... $ 7861.27 
Offprints of Dec. 1952 to Sept.1953 . ...... 1071.07 8932.34 
{ 
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OFFICE OF THE SECRETARY-TREASURER 


Income 
Subscriptions—1952 . $ 286.50 
Back dues and subscriptions .............-+-+8-. 35.00 
Reprints and Back issues... ...... $ 67.60 
Striplist and adjustments. . ....... 34.14 
Less credits taken for overpayments made 
Expenses 
Directory General Total 
Stationery and supplies ... . 33.86 67.42 101.28 
Special services ........ 5.63 75.75 81.38 
Depreciation ......... 230.91 230.91 
Total Expense ..... $1356.45 $6773.31 $ 8129.76 
Excess of Expense over Income ............. $ 1201.02 
ASSETS AND LIABILITIES AS OF DECEMBER 31, 1953 
Assets 
Cash on hand (including $35.00 petty cash). . . 2... $ 2938.91 
1953 dues and subscriptions in transit ........ 452.63 
$ 3622.43 
Liabilities 
Dues and Subscriptions for 1954... 2... 2... $ 1871.00 
Surplus of assets over liabilities . . 1751.43 
$ 3622.43 
Surplus of assets over liabilities, Dec. 31st, 1952 $ 2635.92 
Surplus of assets over liabilities, Dec. 31st, 1953. $ 1751.48 
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WNAR. The Western North American Region met on June 19 at 
the California Institute of Technology in Pasadena. A joint session 
with the Institute of Mathematical Statistics considered social and 
medical applications with papers by Lester Breslow on “Occupations 
and cigarette smoking as factors in lung cancer” and by A. W. Marshall 
on “A study of the onset of mental disease from admission data.’”’ The 
meeting then featured a special invited address by E. C. Hammond on 
“The problem of establishing cause and effect relationships in the etiology 
of chronic diseases.’””’ The afternoon program of contributed papers, by 
D. C. Joseph on ‘Morphological variations in a small marine fish 
(Bairdiella icistius) planted in the Salton Sea” and by John Radovich on 
“Estimating the size of the population of the Pacific sardine (Sardinope 
caerulea)”’, was followed by an informal symposium on methods of teach- 
ing biostatistics to medical students, nurses, sanitarians, zoologists, 
ichthyologists, and other ologists with F. H. Weymouth presiding and 
M. Elveback, C. E. Hopkins, L. E. Moses, W. J. Dixon and D. G. 
Chapman as participants. 

The annual business meeting of the Region was also held on June 19 
with the Regional President, D. G. Chapman, in the chair. After a 
discussion of meetings proposed for December and the next spring, the 
following were elected as Regional Officers for 1955: President—W. J. 
Dixon, Secretary-Treasurer—Elizabeth Vaughan, Regional Committee 
members—D. G. Chapman and B. M. Bennett. 

A meeting of the WNAR of The Biometric Society will be held on 
the Berkeley Campus of the University of California, December 27-29, 
1954, in conjunction with meetings of the American Association for 
the Advancement of Science, the first part of the Third Berkeley 
Symposium, the Institute of Mathematical Statistics, and the American 
Association for the Advancement of Science. Symposiums on the Design 
of Experiments in fisheries; Statistics in Biology and Genetics; and 
Statistics in Medicine and Public Health are being organized. Prof. 
Lincoln Moses of Stanford University, California, is chairman of the 
program committee. 

Région pour la Belgique et le Congo Belge a eu le grand plaisir de 
recevoir le 7 juillet dans les locaux de la Fondation Universitaire, Bru- 
xelles, M. Phillipe F. Bourdeaux, Ph.D., Ingénieur des Eaux et Foréts 
A.I.G.’, Assistant-Professor au North Carolina State College, Raleigh, 
U.S.A. Le sujet du discours de Dr. Bourdeaux etait ‘““Methodes bio- 
metriques dans l’analyse de la vegetation.”” La Conférence etait suivie 
d’une discussion active. ; 
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